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CONVERSION  FACTORS,  U.S.  CUSTOMARY  TO  METRIC  (SI)  UNITS  OF  MEASUREMENT 


U.S.  customary  units  of  measurement  used  in  this  report  can  be  converted  to 
metric  (SI)  units  as  follows: 


Multiply 

by 

To  obtain 

inches 

25.4 

millimeters 

2.54 

centimeters 

square  inches 

6.452 

square  centimeters 

cubic  inches 

16.39 

cubic  centimeters 

feet 

30.48 

centimeters 

0.3048 

meters 

square  feet 

0.0929 

square  meters 

cubic  feet 

0.0283 

cubic  meters 

yards 

0.9144 

meters 

square  yards 

0.836 

square  meters 

cubic  yards 

0.7646 

cubic  meters 

miles 

1.6093 

kilometers 

square  miles 

259.0 

hectares 

knots 

1.852 

kilometers  per  hour 

acres 

0.4047 

hectares 

foot-pounds 

1.3558 

newton  meters 

millibars 

1.0197  x  10"3 

kilograms  per  square  centimeter 

ounces 

28.35 

grams 

pounds 

453.6 

grams 

0.4536 

kilograms 

ton,  long 

1.0160 

metric  tons 

ton,  short 

0.9072 

metric  tons 

degrees  (angle) 

0.01745 

radians 

Fahrenheit  degrees 

5/9 

Celsius  degrees  or  Kelvins1 

*30  obtain  Celsius 

(C) 

temperature  readings 

from  Fahrenheit  (F)  readings, 

use  formula:  C  ■  (5/9) 

(F  -32). 

To  obtain  Kelvin 

(K) 

readings,  use  formula: 

K  -  (5/9)  (F  -32)  +  273.15. 
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A  NUMERICAL  MODEL  TO  SIMULATE  SEDIMENT  TRANSPORT 
IN  THE  VICINITY  OF  COASTAL  STRUCTURES 
by 

Marc  Perl  in  and  Robert  G.  Dean 
I.  INTRODUCTION 


1.  General . 

The  need  for  reliable  predictions  of  shoreline  response  to  man-made  or 
natural  modifications  is  increasing  due  to  environmental  concerns  and  the 
rising  cost  of  remedial  measures.  The  capability  of  numerical  modeling  in 
addressing  problems  of  shoreline  response  has  advanced  with  improvements  in 
wave  climatology,  programs  to  better  understand  sediment  transport 
relationships,  and  improvements  in  numerical  modeling.  In-situ  and  remote 
sensing  technology  for  the  measurement  of  directional  wave  characteristics 
has  been  developed  and  applied,  primarily  within  the  last  two  decades.  In 
addition  to  providing  the  necessary  climatology,  the  resulting  measurements 
have  provided  the  basis  for  evaluation  and  refinement  of  directional  wave 
prediction  procedures.  Studies  such  as  the  Channel  Islands  Harbor  Longshore 
Sand  Transport  Study  (Bruno,  et  al.,  1981)  and  the  Nearshore  Sediment 
Transport  Study  (NSTS)  (Gable,  1979)  have  yielded  a  better  understanding  of 
surf  zone  dynamics  and  the  resulting  sediment  transport.  The  increased 
capacities  of  large  computers  and  reduced  computing  costs  combined  with 
improved  numerical  modeling  algorithms  have  resulted  in  an  extremely 
promising  potential  for  the  numerical  modeling  of  shoreline  problems. 

Although  it  is  doubtful  that  numerical  modeling  will  ever  replace 
completely  the  use  of  movable-bed  physical  models,  the  former  type  offers 
many  advantages.  The  modeling  of  shoreline  response  is  somewhat  analogous  to 
the  problem  of  simulating  storm  surges  in  the  coastal  zone  in  which  the  scale 
effects  and  measurement  difficulties  essentially  preclude  physical  modeling. 
For  shorelines,  the  scale  effects  inherent  in  modeling  sediment  are  well 
recognized  and  the  costs  of  representing  a  substantial  length  of  shoreline 
may  be  prohibitive.  The  laboratory  representation  of  a  realistic  wave 
climate  is  at  the  forefront  of  technology. 

The  investigation  of  shoreline  response  can  best  proceed  by  several 
approaches,  with  each  approach  selected  for  the  particular  strengths  which  it 
offers.  Field  programs  are  costly,  usually  because  of  the  considerable 
equipment  and  the  extensive  time  required,  but  these  programs  are  essential 
for  quantifying  the  values  of  constants  or  parameters,  the  forms  of  which  may 
be  available  from  laboratory  measurements  or  theoretical  considerations. 
Laboratory  studies  occupy  a  special  niche  by  allowing  the  wave  conditions  and 
independent  variables  to  be  controlled  readily,  experiments  to  be  repeated, 
and  selected  measurements  to  be  conducted.  Although,  as  noted  before,  scale 
effects  are  present  in  laboratory  measurements  of  sediment  transport,  the 
physics  governing  the  process  should  be  the  same.  However,  the  relative 
magnitudes  of  suspended  versus  bedload  transport  in  the  laboratory  and  field 
may  differ.  Laboratory  studies  can  also  provide  an  excellent  base  for 
evaluating  certain  aspects  of  a  numerical  model,  including  wave  refraction 
and  diffraction  and  the  resulting  shoreline  patterns  due  to,  for  example,  the 
placement  of  a  littoral  barrier.  Numerical  modeling  offers  the  capability  to 
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incorporate  all  the  tydrodynamlc  wave-surf  zone  and  sediment  transport 
knowledge  that  is  available  from  laboratory  and  field  studies.  Numerical 
modeling  has  the  potential  of  providing  accurate  predictions  of  shoreline 
response  to  various  structural  and  nourishment  alternatives.  Additionally, 
the  possibility  exists  of  employing  numerical  models  and  available  field 
measurements  to  learn  more  about  sediment  transport  mechanisms.  In  this 
latter  mode,  various  candidate  mechanisms  or  coefficients  would  be  evaluated 
by  determining  the  best  match  between  measured  and  predicted  shorelines  and 
the  bathymetry.  Generally,  this  mode  would  require  high-quality  measurements 
of  the  forcing  function  (waves  and  nonwave-related  currents)  and  the 
associated  response  (sediments)  as  well  as  the  knowledge  of  appropriate 
conditions  at  the  boundaries  of  the  model. 

The  present  report  documents  the  development  and  application  of  an 
n-line  numerical  model  to  investigate  bathymetric  response  to  time-varying 
wave  conditions  and  shoreline  modification.  The  model  includes  both 
longshore  and  onshore-offshore  sediment  transport.  Based  on  laboratory 
results,  a  new  distribution  of  longshore  sediment  transport  across  the  surf 
zone  Is  used.  The  wave  climate  is  specified  on  the  model  boundaries  which  do 
not  need  to  extend  to  deep  water.  Efficient  algorithms  are  employed  for 
representing  wave  refraction  and  diffraction.  The  equation  of  sediment 
continuity  and  transport  are  solved  by  a  completely  implicit  algorithm  which 
allows  a  large  time-step.  Specified  sediment  transport  values  or  specified 
contour  positions  can  be  accommodated  at  the  model  boundaries.  The  model  is 
suitable  for  investigating  the  shoreline  response  to  a  variety  of 
modifications  such  as  one  or  more  groins,  terminal  structures,  structures 
with  variable  permeabil ity,  and  beach  nourishment  with  or  without  terminal 
structures. 


Objectives. 


The  objectives  of  the  present  study  include  (a)  the  documentation  of 
state-of-the-art  models,  (b)  the  development  and  documentation  of  an  improved 
model  which  Includes  the  capability  to  represent  n-contour  lines  and  (c)  the 
application  cf  the  model  to  several  relevant  coastal  engineering  problems. 


II.  BACKGROUND 


This  discussion  describes  significant  contributions  which  either  address 
numerical  modeling  of  shorelines  directly  or  provide  Improved  capability  for 
modeling. 

1.  Wave  Refraction  (Noda,  1972). 

Noda  developed  an  algorithm  for  solving  the  following  steady  state 
equation  for  wave  refraction 


?  x  t  -  0  (1) 


in  which?,  the  horizontal  vector  differential  operator,  and  7,  the  wave 
number,  are  defined  in  terms  of  their  components  as 
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*  1  ->  +  1  J 

ax  J  ay 


t 


■t  -  ,  -t  , 

=  i  kx  +  j  ky 


(3) 


where  T  and  j  are  the  unit  vectors  in  the  x  and  y  directions  respectively. 
Equation  (1)  can  be  expressed  as 


a(ksine)  _  a(kcose) 

ax  ”  ay  1  ' 

in  which  e  is  the  direction  of  the  vector  wave  number  relative  to  the  x-axis 
and  k  denotes  | k] .  Noda  expanded  Equation  (4)  to  the  following  form 


k  cose  +  sine  ^  =  -k  sine  ~  +  cose  0  (5) 

9  k  9  k 

Since  —  and  — -  are  known  from  the  angular  frequency  o,  the  water  depth 
dX  dy 

h,  and  the  dispersion  equation 

-  g  k  tanh  kh  (6) 

Equation  (5)  can  be  solved  numerically,  although  there  are  problems  of 
directional  stability.  The  primary  advantage  of  Equation  (5)  is  that  it 
allows  the  wave  direction  e  to  be  determined  on  a  specified  grid,  compared  to 
unspecified  locations  that  would  be  obtained  by,  for  example,  wave  ray 
tracing. 

2.  Crenulate  Bays  (LeBlond,  1972). 

LeBlond  attempted  to  model  the  evaluation  of  an  initially  straight 
shoreline  between  two  headlands  into  a  crenulate  bay.  The  model  constitutes 
a  one-line  (shoreline)  representation.  The  transport  equation  employed 
related  the  total  sediment  transport  to  total  water  transport  in  the  surf 
zone  as  predicted  by  the  formulation  provided  by  Longuet-Higgins  (1970).  The 
initial  shoreline  patterns  resemble  crenulate  bays  in  nature;  however,  the 
predictions  were  found  to  be  unstable  for  reasonably  long  periods  of  computa¬ 
tional  time  and  did  not  approach  a  realfstic  planform. 

3.  Crenulate  Bays  (Rea  and  Komar,  1975). 

Rea  and  Komar  employed  a  rather  Ingenious  system  of  orthogonal  grid 
cells  to  provide  a  cell  which  locally  is  displaced  perpendicular  to  the 
general  shoreline  orientation.  A  one-line  representation  was  employed.  A 
simple  and  approximate  representation  of  wave  diffraction  was  employed. 
Although  the  model  yielded  reasonable  results  for  the  examples  presented,  the 
unique  coordinate  system  would  not  be  suitable  for  a  general  model  as  the 
coordinate  system  must  be  "tailored"  to  some  degree  to  conform  to  the 
expected  shoreline  configurations. 


9 


4.  General  One-line  Shoreline  Model  (Price,  Tomlinson,  and  Willis,  i'iU). 


Price,  Tomlinson,  and  Willis'  formulation  consists  of  the  sediment 
continuity  equation  and  the  total  sedime:  t  transport  equation 


Q 


s 


0.70  (nC)bsinob  cosa^ 


(7) 


in  which  E  represents  the  wave  energy  density,  (nC)  the  group  velocity, 
a  the  angle  between  the  breaking  wave  front  and  the  shoreline,  yu  the 
specific  weight  of  water,  p  the  in-place  sediment  porosity,  and  Ss  the 
specific  gravity  of  the  sediment  relative  to  the  water  in  which  it  is 
immersed.  The  subscript  "b"  represents  values  at  breaking. 

Two  formulations  were  presented  by  Price,  Tom! imson,  and  Willis  (1972). 

In  the  first,  Equation  (7)  was  substituted  into  the  continuity  equation  and 
the  results  cast  into  a  finite-difference  form.  In  the  second,  the  two 
equations  were  employed  separately.  The  latter  formulation  was  selected  due 
to  its  simplicity  and  used  for  the  results  presented. 

Computations  were  carried  out  for  the  case  of  beach  response  due  to  the 
placement  of  a  long  impermeable  barrier.  The  total  sediment  transport 
equation  by  Komar  (1969)  was  used  and  the  planform  was  calculated  at 
successive  times.  Refraction  was  apparently  not  accounted  for  in  the 
numerical  model.  To  verify  the  computations,  a  physical  model  study  was 
carried  out  for  the  same  conditions  using  crushed  coal  as  the  modeling 
material.  The  comparison  was  interpreted  as  good  for  up  to  3  hours;  however, 
for  greater  times,  substantial  differences  occurred  and  these  were  inter¬ 
preted  as  being  due  to  wave  refraction  not  being  represented.  The  crushed 
coal  was  supplied  to  the  model  at  the  updrift  end  at  a  rate  based  on  the 
Komar  equation,  and  the  results  were  interpreted  as  substantiating  this 
relationship.  However,  the  updrift  end  of  the  model  beach  receded  substan¬ 
tially  both  in  the  numerical  and  physical  models.  In  the  physical  model, 
this  can  only  be  interpreted  as  due  to  the  Komar  equation  predictions  being 
less  than  the  actual  transport  rate,  possibly  due  to  the  low  specific  gravity 
(1.35)  of  the  crushed  coal.  The  predicted  recession  of  the  updrift  beach  is 
puzzling,  although  it  could  be  due  to  a  problem  in  properly  representing  the 
updrift  boundary  condition. 

Other  one-line  models  for  shoreline  changes  in  the  vicinity  of  coastal 
structures  were  developed  by  LeMehaute  and  Soldate  (1977)  and  Perlin  (1978). 
Perlin  also  developed  a  two-line  model  formulation,  with  one-line  represent¬ 
ing  the  shoreline  and  the  second  the  offshore.  Dragos  (1981)  developed  an 
n-line  nodel  for  bathymetric  changes  due  to  the  presence  of  a  littoral 
barrier. 
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III.  THE  NUMERICAL  MODEL 


1.  Description. 

There  are  several  methods  of  modeling  bathymetric  changes  due  to  the 
presence  of  a  littoral  barrier.  An  attempt  can  be  made  to  either  model  the 
complete  hydrodynamics  and  the  resulting  sediment  transport  or  model  using  a 
combination  of  analytical  and  empirical  sediment  transport  equations.  The 
second  method  was  chosen  due  to  past  relative  success. 

At  least  two  methods  of  employing  sediment  transport  equations  exist:  a 
fixed  longshore  and  cross-shore  grid  system  where  the  depth  is  allowed  to 
vary  or  a  fixed  longshore  and  depth  system  where  the  cross-shore  distance  is 
allowed  to  change.  Although  it  may  seem  somewhat  awkward,  the  latter  system 
was  chosen  for  the  model.  This  method  allows  the  modeler  to  think  of 
bathymetric  changes  due  to  a  littoral  barrier  in  terms  of  the  effect  on  the 
contours;  i.e.,  the  contour  realinement  due  to  the  structure's  presence  is 
observed.  One  limitation  of  this  approach,  at  least  as  it  was  applied  here, 
is  that  each  depth  contour  must  be  single-valued;  it  is  not  possible  to 
represent  bars. 

The  next  step  in  formulating  the  model  was  choosing  the  specific 
representation  of  the  bathymetry.  The  model  is  an  n-line  representation  of 
the  surf  zone  in  which  the  longshore  direction  x  is  divided  into  equal 
segments  each  ax  in  length.  The  bathymetry  is  represented  by  n-contour 
lines,  each  a  specified  depth,  which  change  in  offshore  location  according  to 
the  equation  of  continuity.  There  are  two  components  of  sediment  transport 
at  each  of  the  contour  lines,  a  longshore  component,  Qx,  and  an  offshore 
component,  Qy.  Figure  lisa  definition  sketch  showing  the  beach  profile 
representation  in  a  series  of  steps  and  the  planform  profile  representation 
and  notations  used. 

Implementation  of  the  sediment  transport  equations  requires  knowledge  of 
the  wave  field  and  the  equilibrium  offshore  profile.  A  discussion  of  the 
refraction  and  diffraction  schemes  follows.  The  equilibrium  profile  is 
introduced  when  it  is  convenient.  As  an  introduction  to  the  logic  used  in 
the  numerical  model,  a  flow  chart  is  presented  in  Figure  2. 

2.  Refraction. 

A  refraction  scheme  compatible  with  variable  Ay’s  was  required  because 
of  the  variable  distance  to  fixed  depth  contours  (as  opposed  to  the  more 
usual  fixed  grid  system  where  a  grid  center  has  a  longshore  and  offshore 
coordinate  with  a  variable  depth).  One  of  the  benefits  of  the  n-line  model 
is  the  ease  with  which  the  response  of  the  contours  to  a  particular  wave  and 
structure  condition  can  be  visualized.  A  fixed  grid  system  and  an 
interpolation  scheme  could  have  been  used  to  obtain  the  wave  field;  however, 
this  would  have  reduced  accuracy  and  increased  computation  time.  The  scheme 
developed  also  saves  computation  time  because  it  does  not  use  differential 
products  terms. 
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(b)  Beach  Plonform  Representation 

Figure  1.  Definition  sketch. 


The  first  of  the  governing  equations  used  is  the  conservation  of  waves 
equation 


,8) 

->  -► 

where  vh  is  the^ horizontal  differential  operator  equal  to  i(a/ax)  +  j ( a/ay ) 
in  which  i  and  j  are  the  unit  vectors  in  the  x  and  y  directions,  respec¬ 
tively,  and  x  is  the  longshore  direction,  with  positive  to  the  right  when 
facing  the  water,  y  the  offshore  direction,  with  positive  seaward,  and  z  the 
vertical  coordinate,  with  positive  defined  as  upwards.  For  the  steady-state 
case,  equation  (8)  yields 


h  'V  -  h  ■ 0 


(9) 


where  kx  and  ky  are  the  wave  number  projections  in  the  respective  directions. 
Defining  as  the  angle  k  makes  with  the  y-axis  positive  in  the  counter¬ 
clockwise  direction,  the  equation  can  be  written  in  final  form  as 


fx  (k  C0S  6)  =  iy  {k  S1’n  9) 


(10) 


where  e  =  a  +  ir  (in  radians).  Noda  (1972)  and  others  have  developed 
numerical  solutions  to  expanded  forms  of  equation  (10).  In  the  present 
study,  equation  (10)  was  initially  central -differenced  in  the  x-direction  and 
forward-di fferenced  in  the  y-direction  with  Snell's  law  used  to  specify  the 
boundary  conditions  on  the  offshore  boundary  and  one  of  the  sides  (i.e.,  the 
side  of  the  wave  angle  approach).  However,  a  numerical  problem  arose.  The 
argument  of  the  arcsine  exceeded  +1.0  for  large  ay/ ax.  To  overcome  this 
problem,  a  dissipative  interface  was  used  on  the  forward-difference  term 
(after  Abbott,  1979).  The  final  finite-di fferenced  form  of  equation  (10)  is 
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where  t  has  been  taken  as  0.25.  The  past  ©"j  and  the  present  0,  j  wave 
angles  are  numerically  averaged  to  give  the  0-j  Newton's  method  is  used 
to  compute  the  wave  number  via  the  linear  wave  theory  dispersion  relation. 

In  addition,  numerical  smoothing  is  used  at  the  conclusion  of  the  wave  field 
calculation.  This  approximates  in  an  ad  hoc  manner  diffractive  effects 
(lateral  transfer  of  wave  energy  along  the  wave)  which  exist  in  nature  but 
have  been  omitted  due  to  use  of  the  equation  for  refraction  (equation  8). 

The  smoothing  routine  is 


1 ,  j 
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The  second  governing  equation  used  in  the  refraction  scheme  is 
conservation  of  energy.  Neglecting  dissipation  of  energy  due  to  friction, 
percolation,  and  turbulence,  this  equation  can  be  expressed  as 


v  -(E  £q)  =  0 


(13) 


where  E  is  the  average  energy  per  unit  surface  area  and  Cg  the  group  velocity 
of  the  wave  train.  Performing  the  operation  indicated  and  replacing  (L  by 
its  components  (Crsin  e)  and  (Crcos  e)  results  in  the  following:  “ 


~  (E  Cr  sin  «)  +  j-  (E  C  cos  e)  =  0 
3x  b  ay 


(14) 


Assuming  linear  theory, 

2 

pgfT 


(15a) 


where  p  is  the  mass  density  of  water,  g  th^agravitational  constant,  and  H 
the  wave  height.  Dividing  the  equation  by  jp,  finite-differencing  and 
weighting  the  forward-di fferenced  term  as  before,  and  solving  for  the  wave 
height,  results  in  the  following: 
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This  equation  is  also  solved  by  iterative  techniques  and  the  and  Hn.  . 
are  averaged  at  the  conclusion  of  each  iteration. 

Cg  is  determined  by  the  linear  wave  theory  relationship 


n  _  C  / i  ,  2kh  i 

CG  ~  2  1  +  (16) 


where  h  is  the  water  depth,  k  the  wave  number,  and  C  the  wave  celerity.  Wave 
height  boundary  conditions  are  input  along  the  same  boundaries  as  the  wave 
angles  using  linear  theory  shoaling  and  refraction  coefficients.  The  e's 
have  been  previously  determined.  In  both  equations  (11)  and  (15)  for  a 
variable  grid  system,  the  points  (i  +  1,  j)  and  (i  -  1,  j)  need  to  be 
determined  (i.e.,  because  the  y  coordinates  are  not  ffxed,  adjacent  values 
with  the  same  subscripts  can  be  farther  or  closer  to  shore,  therefore 
Interpolation  must  be  used).  The  actual  values  are  found  by  searching  the 
(i  +  1)  and  (i  -  1)  cross-shore  lines,  finding  the  adjacent  values  in  the 
positive  and  negative  y-direction,  and  interpolating  to  determine  the  value. 


3.  Diffraction. 

The  diffraction  solution  (in  the  lee  of  the  structure)  used  in  the  model 
is  based  on  the  method  of  Penny  and  Price  (1952).  Assumptions  used  in  this 
method  include  a  semi-infinite  breakwater,  which  is  infinitesimally  thin, 
linear  wave  theory  and  constant  depth.  A  definition  sketch  for  wave 
diffraction  is  shown  in  Figure  3.  The  quantity  THETAO  represents  the  angle 
of  wave  incidence  relative  to  the  jetty  axis,  ANGLE  represents  the  angle  from 
the  jetty  at  the  point  where  the  diffraction  coefficient  is  to  be  computed, 
and  RAD  is  the  radial  distance.  The  radial  distance  is  then  cast  into  a 
dimensionless  parameter,  RHOND  (=  2*  RAD/L),  where  L  is  the  wavelength.  This 
is  equivalent  to  multiplying  the  radial  distance  by  the  wave  number  k. 

The  diffraction  coefficient  AMP  is  expressed  as  the  modulus  of  the 
diffracted  wave 

AMP  =  (Sum  l)2  +  (Sum  2)2  (17) 

where 

Sum  1  =  [cos  (RHOND  (cos  (ANGLE-THETAO)))  .  (g  (1.0  +  Cp  +  S))]  + 

[sin  (RHOND  (cos  (ANGLE-THETAO)))  .  (-  \  (S  -  Cp))]  + 

[cos  (RHOND  (cos  (ANGLE+THETAO) ) )  .  (\  (1.0  +  Cp  +  S))]  + 

[sin  (RHOND  (cos  (ANGLE+THETAO)))  .  (j  -  (S  -  Cp))]  (18) 


Sum  2  =  [cos  (RHOND  (cos  (ANGLE-THETAO)))  .  (-  \  (S  -  Cp))]  + 

[sin  (RHOND  (cos  (ANGLE-THETAO)))  .  (]>  (1.0  +  Cp  +  S))]  + 

[cos  (RHOND  (cos  (ANGLE+THETAO)))  .  (-  \  (S  -  Cp) ) ]  + 

[Sin  (RHOND  (cos  (ANGLE+THETAO)))  .  (\  (1.0  +  Cp  +  S))]  (19) 


In  Equations  (18)  and  (19),  Cp  and  S  represent  Fresnel  integrals  and  are 
computed  in  the  model  by  means  of  an  approximation  after  Abramowitz  and 
Stegun  (1965). 

Having  obtained  AMP,  the  wave  height  at  the  location  in  question  is 
simply  the  product  of  the  specified  partially  refracted  incident  wave  height 
and  AMP.  The  angle  of  the  wave  crest  is  computed  assuming  a  circular  wave 
front  along  any  radial;  this  angle  is  then  refracted  using  Snell's  law. 

Throughout  the  refraction  and  diffraction  schemes,  the  local  wave 
heights  are  limited  by  the  value,  0.78  x  depth. 
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Figure  3.  Definition  sketch  for  wave  diffraction. 
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4.  Sand  Transport  Model. 

a.  Governing  Equations.  Three  basic  equations  are  used  to  simulate  the 
sediment  transport  and  bathymetry  changes  according  to  the  wave  field.  The 
equation  of  continuity 

aqx  aq 

+  +  — *  =  0  (20) 
ax  ay 


requires  as  input,  knowledge  of  the  longshore  and  cross-shore  components  of 
sediment  transport.  The  total  transport  alongshore  has  been  measured  by 
several  investigators  and  many  equations  exist;  however,  the  distribution  of 
the  transport  across  the  surf  zone  is  not  well  known.  Fulford  (1982)  based 
on  laboratory  data  from  Savage  (1959),  developed  a  distribution  of  longshore 
sediment  transport  across  the  surf  zone  for  the  case  of  straight  and  parallel 
contours.  Fulford's  use  of  Savages  experiment  was  based  on  two  assumptions: 
1)  the  structure  must  be  a  total  littoral  barrier  and  2)  onshore-offshore 
sediment  transport  could  be  neglected.  Test  5-57  was  chosen  because  the  two 
criteria  were  nearly  met.  Savage  reported  that  the  groin  acted  as  a  total 
littoral  barrier  for  the  first  35  hours  of  the  test  (i.e.,  no  bypassing 
occurred  prior  to  35  hours).  This  does  not  mean  that  no  onshore-offshore 
transport  occurred  because  as  the  profile  steepens  on  the  updrift  side, 
onshore-offshore  transport  does  occur.  However,  it  was  assumed  to  be 
negligible.  In  addition,  the  initial  profile  had  been  molded  to  an 
equilibrium  profile  via  150  hours  of  waves.  Thus,  the  two  criteria  required 
to  develop  an  inferred  longshore  distribution  of  sediment  transport  were 
nearly  satisfied.  This  distribution  is  shown  as  a  dashline  in  Figure  4.  The 
smaller  "maximum"  is  believed  to  be  an  extraneous  effect  of  a  groin  downdrift 
from  the  location  in  the  experiment  where  the  data  were  taken.  Therefore, 
this  feature  was  replaced  by  a  monotonical ly  decreasing,  smooth  curve  as 
shown  by  the  "altered"  curve.  To  analytically  represent  this  distribution,  a 
function  of  the  following  form  was  chosen 

qx  (y)  «  (B)  (y ) n~ 1  e'{y)"  (21) 

This  type  of  equation  is  convenient  because  it  is  easily  integrable,  and  by 
properly  choosing  the  constant,  B,  the  integral  of  the  equation  from  zero  to 
infinity  can  be  required  to  equal  a  particular  value.  This  too  is  highly 
desirable  because,  as  was  done  in  the  model,  the  integral  is  set  equal  to  one 
and  then  multiplying  by  the  value  of  the  well-known  longshore  transport 
equation,  the  value  of  the  transport  at  any  location  across  the  surf  zone  can 
be  determined.  Further  investigation  suggested  a  value  of  n  =  3  to  produce  a 
curve  similar  to  Fulford's  curve.  A  more  general  form  of  the  equation  which 
allows  more  flexibility  and  curve  fitting  is 


qx(y)  =  B(y  +  a)2 


(22) 
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Figure  4.  Distribution  of  sediment  transport  across  the  surf  zone. 


where  yb  =  distance  to  the  point  of  breaking 

a  =  constant  to  allow  sediment  transport  above  mean  water  line  (MWL) 
(swash  transport  or  transport  in  region  of  wave  setup)  to  be 
represented 

c  =  a  constant  establishing  the  width  of  the  curve  (to  be  determined) 


B  = 


"Tl 

c  ys 


[causes 


/ 


Qx  «'y)  dy  =  1.0) 


Based  on  Fulford’s  (1982)  results  and  considering  a  to  be  proportional 
to  the  breaking  height  divided  by  the  beach  slope,  the  constant  of  propor¬ 
tionality  was  determined  to  be  unity;  i.e.,  a  =  hb/(ah/ay).  Using  equation 
(22)  and  a  digitized  version  of  the  curve  shown  in  Figure  4,  a  nonlinear  least 
squares  regression  was  carried  out  to  determine  the  value  of  c.  A  Taylor's 
series  expansion  of  the  form 


fk  +  1  (c,y)  =  fk(c,y)  +  £ 


AC 


(23) 


where  k  and  k  +  1  represent  the  number  of  the  iteration  carried  out.  Least 
squares  regression  minimizes  the  square  of  the  difference  between  observed 
and  predicted  values  with  respect  to  a  change  in  the  parameter  being 
computed,  or 


a 

aUc] 


U  t 


OBS 


-  (fK(c,y)  + 


(24) 


where  foBS  represents  the  observed  values,  which  in  this  case  is 
qx(y)oBS-  Carrying  out  the  differentiation  indicated  and  manipulating 
terms,  ac  can  be  solved  in  terms  of  known  quantities. 

An  iterative  procedure  was  then  used  by  updating  the  values  of 
fk(c,y),  af/ac,  and  c  until  an  acceptably  small  change  in  c  results.  For 
the  data  herein,  the  value  of  c  was  determined  to  be  1.25.  The  final  form  of 
sediment  transport  of  a  y  location  in  the  surf  zone  results  for  a  shoreline 
with  straight  and  parallel  contours,  as 


qx(y)  - 


(yta)2  ,-[(>*«)/(!. 25  y„)r 


(25) 
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This  equation,  which  is  also  presented  in  Figure  4,  predicts  the  relative 
transport  at  point  y.  To  obtain  the  fraction  of  transport  between  two  y 
coordinates,  the  integral  of  equation  (25),  from  yj  to  yg,  must  be 
used. 


qx(y)4y  -  e  -[(>1  *  al/(1-25 

_e  -Hy2+  a )/{ 1 . 25  yb)]3 


(26) 


QX[ND]  is  dimensionless;  therefore,  to  compute  a  value  in,  say,  cubic  feet  per 
second,  it  must  be  multipled  by  the  total  transport  along  a  perpendicular  to 
the  shoreline  obtained  from  the  total  longshore  transport  equation  used  in 
the  model 


Q  =  C'  Hb5/2  sin  (2  ab) 


(27) 


See  Appendix  A  for  a  discussion  of  the  constant  C'.  It  is  noted  that  the 
transformation  of  qx(y)  to  qx(h)  can  be  effected  by  multiplying  by  the 
one-dimensional  Jacobian  (ay/a h).  This  latter  form  (qx(h))  is  more  useful 
here  because  the  present  model  simulates  the  changes  in  contour  position  (ay) 
rather  than  changes  by  depth  (ah). 

In  the  numerical  model,  Qx  (I,J)  (see  Fig.  1)  is  determined  using 
equation  (26)  except  for  the  shoreline  contour,  J=l,  and  the  farthest 
offshore  contour  simulated,  J  =  JMAX.  The  shoreline  contour  longshore 
transport,  Qx  (1,1),  in  order  to  include  swash  transport,  uses  equation 
(16);  however,  the  first  term  is  set  equal  to  1.0.  The  seawardmost  contour 
transport,  Qx  (I, JMAX),  in  order  to  include  any  longshore  transport  not  yet 
accounted  for,  neglects  the  second  term  of  equation  (26)  (i.e.,  it  accounts 
for  transport  from  y(I,JMAX)  to  infinity).  The  dimensionless  numbers  are 
then  multiplied  by  Q  determined  from  equation  (27).  This  method  is  based  on 
parallel  contours  which  may  not  exist.  In  order  to  compensate  for  the 
nonparallel  nature  of  the  contours  (note  that  refraction  does  account  for  it 
as  far  as  the  wave  field  is  concerned),  the  term  sin  (2a|>)  of  equation  (27) 
is  replaced  by  sin  (2ai)  shoreward  of  the  breakpoint,  where  a[_  represents 
the  angle  between  the  "local"  wave  angle  and  the  "local"  contour.  It  can  be 
argued  that  for  a  spilling  breaker,  the  remaining  surf  zone  at  any  point 
"sees"  a  total  transport  similar  to  equation  (27),  where  05  and  Hj,  are 
the  local  values.  The  problem  is  that  the  constant  of  proportionality  was 
determined  for  the  entire  surf  zone  and  for  nearly  straight  and  parallel 
contours.  This  not  being  the  case,  the  equation  was  altered  on  intuitive 
grounds  to  reflect  the  fact  that  the  contours  are  no  longer  straight  and 
parallel . 
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The  second  input  required  by  the  continuity  equation  to  predict  the 
bathymetric  changes  is  the  cross-shore  sediment  transport.  The  governing 
equation  for  onshore-offshore  transport  (after  Bakker,  1968)  is 


AX 


r .  .  i 


(28) 


_5 

where  C^r-  is  an  activity  factor  (inside. the  surf  zone  =  10  feet  per  second 
for  theuprototype  simulation  herein,  1(j  feet  per  second  for  the  physical 
model  simulation)  (see  App.  A.  for  a  discussion)  and  W^gfiJ)  is  the 
positive  equilibrium  profile  distance  between  y(i,j)  ana  y(i,j-l).  determined 
from  the  equilibrium  profile  used  in  the  numerical  model  h  =  Ay^/3  (Dean, 
1977).  See  Appendix  A  for  discussion  of  the  value  of  A.  The  physical 
interpretation  of  equation  (28)  is  that  as  this  profile  steepens  (flattens), 
sediment  is  transported  offshore  (onshore). 


b.  Methods  of  Solution.  Three  separate  finite-difference  techniques 
were  used  to  solve  the  equations: 

(1)  Explicit  longshore-continuity  and  explicit  cross-shore 
continuity; 

(2)  Implicit  longshore-continuity  and  explicit  cross-shore 
continuity  for  half  a  time-step  then  vice  versa;  and 

(3)  Implicit  longshore-cross-shore  continuity. 

An  explicit  formulation  was  first  developed  which  used  the  refraction 
scheme,  the  distribution  of  longshore  sediment  transport  across  the  surf 
zone,  and  the  onshore-offshore  sediment  transport  equation.  Problems  in 
addition  to  the  usual  ones  which  are  encountered  with  explicit  methods  (e.g., 
computation  time  and  cost)  were  immediately  realized.  In  the  explicit 
method,  both  transport  computations  are  based  on  the  former  values  of  the 
contour  locations  and  are  completely  uncoupled.  Stability  of  an  explicit 
scheme  requires  a  small  time-step.  In  addition,  the  noncoupled  nature  of  the 
equations,  in  some  cases,  resulted  in  crossing  of  the  contours  due  to  the 
transport  computed. 

It  is  logical  to  assume  that  an  implicit  formulation  of  the  longshore 
transport  equation  used  as  input  to  the  continuity  equation  along  with  the 
explicit  onshore-offshore  transport  component  would  help  the  numerical 
stability  (on  the  other  half  time-step,  the  longshore  component  would  be 
computed  explicitly  and  the  onshore-offshore  transport  equation  would  be 
solved  implicitly  with  the  continuity  equation).  Although  this  scheme  would 
be  superior  to  the  explicit  procedure,  it  still  would  be  susceptible  to 
crossing  contours.  It  should  be  noted  that  the  magnitude  of  the  coefficient 
used  in  the  onshore-offshore  equation  is  very  important  to  the  extent  tha* 
the  simulation  models  natural  phenomena.  If  the  coefficient  is  very  small  or 
vanishes,  sediment  will  not  move  offshore  and  contours  will  cross  because  of 
the  variation  in  the  distribution  of  longshore  sediment  transport  across  the 
surf  zone.  If  the  coefficient  is  too  large,  the  onshore-offshore  transport, 
may  become  large  enough  that  on  a  particular  time  step,  an  offshore  contour 


\ 
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would  move  too  far  shoreward,  thereby  crossing  an  inshore  contour  or  vice 
versa.  Once  the  contours  cross,  not  only  does  the  bathymetry  become 
unrealistic,  but  mathematically,  the  equation  which  computes  the  longshore 
distribution  across  the  surf  zone  changes  signs  at  some  locations  and  the 
entire  model  becomes  physically  unrealistic. 

To  circumvent  these  problems,  an  implicit  scheme  that  simultaneously 
solves  the  three  governing  equations,  was  developed.  Utilizing  equation  (26), 
and  the  one-dimensional  Jacobian  ( Ay/Ah)  to  convert  to  Qx(h),  the  total 
longshore  transport  equation  (27),  the  following  equation  is  obtained. 


Qx(i,j)  represents  the  sediment  transport  between  depths  h(i,j)  and  h(i,j-l) 
(see  Fig.  1).  The  term  in  brackets  represents  the  normalized  distribution  of 
longshore  transport  between  h(i,j)  and  h(i,j-l);  e  is  the  averaged  wave  angle 
at  the  location  of  Qx(i,j)  and  ac  is  the  local  contour  orientation  angle. 
Defining  everything  except  sin  (2e  -  2ac)  as  v(i,j)  and  using  a  superscript 
to  denote  a  time  step,  this  equation  can  be  written 


v.  .  sin  (2e  -  2a"+1) 

I  j  J  V* 


(30) 


The  assumption  has  been  made  that  the  wave  field  (H  and  e)  do  not  vary  during 
the  bathymetric  changes  over  the  time-step.  Using  the  following  trigonometric 
identities, 


sin 

(2a  - 

2b)  =  sin  2a  cos  2b  -  cos  2a  sin  2b 

(31a) 

cos 

2a 

=  2  cos2  a  -  1 

(31b) 

sin 

2a 

=  2  sin  a  cos  a 

(31c) 

and  recognizing  that  the  following  expression  is  an  approximation 


1  ,..n+l 


n+1 


sfn 


-  .  v"TA  +  y"  .  y"  } 

-  I  lyi,j  yi-l,J  yiJ  ^-lj* 


(32) 
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along  with  assuming  that  the  change  in  the  denominator  is  small  for  a 
reasonable  time-step  (the  numerator  has  been  averaged  over  the  nth  and  n  + 
1th  time-steps),  equation  (30)  results  in 


+  (S3) 


1  J 


n+1 

1-1 J 


=  (RHSl)J  . 

•  J  J 


(33) 


where 


(S3).  .  =  (*)  (v.  .)  cos  (2e)  (2  cos  a) 
i.j  2  i,j  c 

{ RHS1 ) ?  *  =  (v.  •)  (2  sin  e  cos  e)(cos2a  - 
i.j  i,j  c 


1 _ 

(**2  *  J  )  1/2 
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Here  it  has  also  been  assumed  that  cos  ac  does  not  change  over  the  time 
step.  Equation  (33)  is  the  final  form  of  the  longshore  sediment  transport 
equation  prior  to  its  use  in  conjunction  with  the  other  equations. 

Averaging  y  values  on  the  n^  and  (n+l)th  time-steps,  equation  (29) 
car.  be  rewritten  as 


=  Const6-  . 

■  *  J 


+  y 


i,j-l 


(34) 


where  Const6(i,j)  =  CoFp(1>j)«  ax.  This  is  the  final  form  on  the 
onshore-offshore  sediment  transport  equation. 


The  equation  of  continuity,  finite-differenced  for  the  nth  and 
(n+l)th  time-steps,  can  be  written  as 


yn-+1-yn. 


At 


J  = 


^  j  q;+1  +q;  -q;+1  <  -q"+1  -q;  ! 

(  i,J  i+l,j  i+l,j  yi,j  yi,j  yi ,  j+1  yi  ,j+li 

(35)  1 


Defining  Rfj  as  l/(2axah),  inserting  equations  (33)  and  (34)  into  equation 
(35),  and  transferring  all  known  quantities  for  the  nth  time-step  to  the 
right-hand  side  of  the  equation  result  in 


n+1 


yi,j  +  (AtRi,j)S3i,jyi,j  "  UtRi,j)S31,jyi-l,j  '  (AtRl,j)S3l+lJy"+lj 

(AtRi,j)S3i+l,jy1,j  '  (AtRi,jConst6i,j)  (2  [  yi j-l  "  yi J  ]) 
(AtRi,jConst6i  J+i)  (1  [  yjj  -  y" j+1  ])  -  (AWARE) (36) 
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Equation  (36)  can  be  rewritten  as 


(1  +  u  .  V  *  Z1  ♦  Z2)  y£»  -  (V)y^}fJ 

'  -  ,22<]*1  '  (A“ARE)i,j 


(37) 


where 


u  ■  AtR,-,jS3f,j 
V  ■  lB1Js3W,i 
Z1  -  (ii)  Ri  ,Const6,  , 

Z2  ■  (“)RlJConst6j>Jtl. 


Equation  (37)  is  a  weighted,  centered  scheme  in  which  y"  j  is  computed 
using  a  weighting  of  itself  and  its  four  adjacent  grid  "neighbors".  The 
weighting  factors  (U,  V,  Zl,  and  Z2)  are  functions  of  the  wave  climate,  the 
slope  between  contours,  and  the  variables  included  in  the  original  formula¬ 
tion.  An  investigation  of  a  small  gridded  system  demonstrated  that  by  writing 
simultaneous  equations,  one  for  each  yi,j,  a  banded  matrix  results.  This 
matrix  can  be  solved  by  LEQT1B,  one  of  tne  available  routines  from  the 
International  Math  and  Statistics  Library  (IMSL).  A  schematic  representation 
of  the  matrix  A  which  results  from  the  matrix  equation  [A][y]  =  [B]  is 
presented  in  Figure  5.  In  this  schematic,  the  large  zeros  represent 
triangular  corner  sections  of  all  zeros  and  the  0...0  represents  bands  of 
zeros,  the  number  of  which  is  dependent  on  the  number  of  contours  simulated 
(the  number  of  zero  bands  between  either  remote  nonzero  bands  and  the 
tridiagonal  nonzero  bands  equals  two  less  than  the  number  of  contours  modeled 
(in  both  the  upper  and  lower  codiagonals  of  the  matrix)).  An  inspection  of 
the  subscripts  in  equation  (29)  yields  the  reason  the  zero  bands  are  required. 
The  more  j  values  (contours)  used,  the  more  y  grids  there  are  along  any 
perpendicular  to  shore.  This  causes  zeros  to  appear  in  the  matrix  between 
bands  as  the  weighting  factors  await  being  used  to  operate  on  yn+1(i-l,j) 
and  yn+^(1+l,j).  For  this  reason,  the  expense  of  simulating  an  increasing 
number  of  contours  is  exponential.  The  LEQT1B  routine,  utilizes  banded 
storage  and  saves  both  storage  and  computation  time;  however,  the  routine  has 
no  special  way  of  handling  the  interior  zero  bands.  One  refinement  which 
would  save  computation  time  would  be  to  develop  an  algorithm  to  solve  and 
store  the  matrix  by  taking  advantage  of  these  inner  zero  bands;  however,  it 
is  beyond  the  scope  of  this  project. 

Of  course,  the  matrix  requires  boundary  values  on  longshore  extrem¬ 
ities  and  on  both  onshore  and  offshore  boundaries.  The  longshore  boundary 
conditions  are  treated  by  modeling  a  sufficient  stretch  of  shoreline  so  that 
effects  of  a  structure's  presence  are  minimal.  The  y  values  along  these 
boundaries  can  therefore  be  fixed  at  their  initial  locations.  In  the 
onshore-offshore  direction,  boundaries  are  treated  quite  differently.  The 
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Note:  Size  of  matrix  full  storaqe  mode 
[ ( IMAX-2) ( J MA X )  x  ( IMAX-2 ) (JMAX ) ] 

Size  of  matrix  banded  storaqe  mode 
[(IMAX-2) (JMAX)  x  (2JMAX  +  1  )] 


Figure  5.  Schematic  representation  of  banded  matrix  if 
not  stored  in  banded  storaqe  mode. 


1 


berm  and  beach  face  are  assumed  to  move  in  conjunction  with  the  shoreline 
position.  The  required  sediment  transport  is  then  computed  by  the  change  in 
position  of  the  shoreline.  The  two  equations  are 


r  n+1 

[yi,l 


(38a) 


(38b) 


The  offshore  boundary  is  treated  by  keeping  yn+^(i,jmax)  (the  contour 
beyond  the  last  simulated  contour)  fixed,  until  the  angle  of  repose  is  exceeded. 
Then,  the  yn+1 ( i , jmax+1 )  is  reset  (at  the  conclusion  of  the  n+1  time-step) 
to  a  position  such  that  the  slope  equals  the  angle  of  repose.  Note  that 
yn+^(i,0)  is  represented  in  the  program  by  YZERO-j . 


There  are  also  no-flow  boundary  conditions  required  at  each  of  the 
structures  being  modeled.  These  are  imposed  on  the  adjacent  y-grid  points  which 
are  located  downdrift  (i.e.,  in  the  shadow  zone)  of  the  structure  and  shoreward 
of  the  structures'  seaward  extremities.  They  are  imposed  by  setting  S3j  j  of 
equation  (33)  and  DISTRij  (the  term  in  square  brackets  in  equation  (29) 
equal  to  zero,  thereby  causing  Qx(i,j)  to  be  zero  (i.e.,  the  no-sediment  flow 
condition).  This  boundary  condition  is  imposed  automatically  for  every 
shore-perpendicul ar  structure. 


It  was  found  that  even  with  the  implicit  formulation,  high  frequency 
oscillations  occurred  in  the  y  values  immediately  updrift  and  downdrift  of  the 
structure.  The  solution  did  not  "blow  up";  however,  on  larger  time-steps 
"sloshing"  (oscillating)  did  occur.  Part  of  this  problem  was  due  to  the 
boundary  condition  at  the  structure  which  had  been  such  that  either  no  sand  was 
allowed  along  a  contour  line  or  the  sand  determined  by  the  equations  was  allowed 
to  be  transported.  Because  of  the  very  large  angle  which  existed  around  the  tip 
of  the  structure  when  a  contour  first  exceeded  the  length  of  the  structure,  very 
large  amounts  of  sediment  transport  were  predicted.  In  the  nature  where  analog 
sand  transport  rather  than  digitized  transport  occurs,  this  does  not  happen. 
Therefore,  the  boundary  condition  was  altered  to  constantly  allow  sand  transport 
around  the  end  of  the  structure  in  proportion  to  that  part  of  the  contour 
representation  which  exceeded  the  structure  (i.e.,  the  transport  was  calculated 
for  the  location  at  tip  of  the  structure  as  if  the  structure  was  not  there  and 
then  a  proportion  of  this  value  was  allowed  to  bypass).  Although  the  transport 
around  the  tip  of  the  structure  is  based  on  the  values  from  the  past  time- step, 
it  more  closely  simulated  the  natural  phenomenon. 


Additionally,  a  dissipative  interface  is  used  on  the  y  values  as  follows: 


i,J 


=  (t)  y 


1-l.j 


+  (1  -  2t)  y.  .  +  (t)  y 
'  >  J 


i+l.j 


(39) 


where  t  was  again  taken  as  0.25.  It  is  noted  that  only  high  frequency 
oscillations  in  y  are  affected  by  the  use  of  equation  (39);  the  total  sum 
of  y  values  is  not  affected.  Also,  in  all  the  dissipative  interface 


schemes  used,  if  a  boundary  point  is  being  computed,  either  a  forward- 
difference  or  a  backward-difference  of  equation  (39)  is  used  (after 
Abbott,  1979): 


(40a) 
(40b) 

IV.  SIMULATIONS  AND  VERIFICATION 


Backward:  y(J  -  My,.u  *  »  -  .  >y,  j 
Forward:  y,  j  =  hly|(y  ♦  0  -  >  >»,  j 


Several  simulations  were  run;  two  were  attempts  at  verifying  the 
numerical  model,  the  others  were  run  to  gain  insight.  Because  a  complete 
data  set  does  not  exist,  only  the  available  data  are  compared.  The  first 
modeling  effort  was  to  simulate  the  physical  model  tests  of  Savage  (1959).  A 
second  set  of  cases  was  run  for  shore-perpendicular  structures.  Next,  an 
effort  was  made  to  model  sediment  transport  in  the  vicinity  of  a  hypothetical 
dredge  disposal  site  in  the  11-  to  14-foot  depths  off  Oregon  Inlet.  Finally, 
the  Channel  Islands  Harbor  Longshore  Transport  Study  (Bruno,  et  al . ,  1981) 
was  modeled.  Bathymetric  changes  were  closely  monitored  during  this  study; 
however,  the  wave  climate  (H,  e,  T)  used  was  determined  from  the  Littoral 
Environmental  Observation  (LEO)  data  and  uncertainties  exist  as  to  the 
accuracy  of  the  data. 

1.  Simulation  of  Savage's  Physical  Model  Tests. 

The  numerical  model  was  used  to  simulate  one  of  the  physical  model  tests 
of  Savage  (1959).  Test  5-57  was  simulated  numerically  for  a  10-hour  period. 
In  this  physical  model,  the  mean  sediment  size  was  0.22  millimeters,  the  wave 
height  averaged  0.25  feet,  the  wave  period  was  1.5  second,  the  wave  angle  was 
30°  (at  a  depth  of  2.3  feet),  and  the  groin  was  approximately  9.5  feet  from 
still  water  to  its  seaward  limit.  Coff  was  he1d  constant  at  10-4  feet 
per  second  throughout  the  profile  for  this  simulation.  The  offshore  profile 
is  presented  in  Savage  (1959).  Figure  6  represents  three  of  the  eight 
contours  simulated.  Note  that  the  initial  0.3-  and  0.5- foot-depth  contours, 
in  the  numerical  representation  are  too  far  seaward  by  approximately  2  feet. 
This  is  due  to  the  h  =  Ay^'4  equation  as  compared  to  the  equilibrium 
physical  model  profile.  Realizing  this,  it  is  the  shape  of  the  contour  which 
must  be  used  as  an  indication  of  the  numerical  model  predictions.  The 
general  trend  of  the  contours  is  similar,  although  the  numerical  model 
contours  are  displaced  farther  seaward  as  expected.  The  major  differences 
are  in  the  diffraction  zone. 

2.  Several  Runs  Using  Shore  Perpendicular  Structures  to  Demonstrate  Effects 
of  Altering  Some  of  the  Pertinent  Parameters. 

In  the  following  simulations,  the  models  were  run  until  their 
near-equilibrium  values  were  achieved.  Coefficient  Coff  w4s  not  a  function 
of  depth  (beyond  the  surf  zone)  but  was  held  constant  throughout  the 
simulated  area.  Important  variables  are  as  shown  in  the  figures.  Only  one 
wave  condition  (H0  *  3  feet,  T  =  7  seconds,  and  a  deepwater  wave  angle  a0 
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Longshore  Grid  Locations  (DX  =  a . 44  ft) 

Figure  6.  Simulation  of  the  physical  model  of  Savage  (1959). 


of  60°)  was  used  as  input  for  all  four  cases.  Case  4.2a  used  an  equilibrium 
shape  factor  A  of  0.0899  and  one  groin.  Case  4.2b  was  similar  to  4.2a  with 
the  only  modification  being,  that  the  A  value  was  changed  to  0.1486.  In  this 
way,  a  direct  comparison  was  made  based  only  on  the  shape  of  the  equilibrium 
profile.  Cases  4.2c  and  4. 2d  used  A-values  of  0.0899  and  0.1486, 
respectively,  but  this  time  three  shore-perpendicular,  evenly  spaced 
structures  were  simulated. 

a.  Comparison  of  Cases  4.2a  and  4.2b.  The  most  obvious  difference  between 
Figures  7  and  8  is  the  volume  of  sand  impounded  updrift  and  eroded  downdrift. 

This  is  due  to  blockage  of  more  of  the  active  transport  zone  in  the  second 

case  (i.e.,  a  shorter  groin  is  required  for  an  equivalent  performance  on  a 
steeper  beach).  The  next  obvious  difference  is  the  size  of  the  perturbation 
which  exists  in  the  offshore  contours.  Clearly,  case  4.2b  is  more  perturbed 
and  this  is  expected  because  larger  offshore  transports  occur  due  to  the 
steepening  on  the  updrift  side.  Conversely,  this  means  less  sediment  is 
initially  bypassed  (and  along  with  the  downdrift  requirement  for  larger 
volumes  of  sand)  causes  larger  erosional  features  in  case  4.2b.  Another 
interesting  feature  is  the  downdrift  fillet  which  occurs  in  the  third, 
fourth,  and  fifth  contours.  The  fillet  is  due  to  the  shape  of  the  sixth 
contour  which  occurs  because  of  the  inability  of  the  wave  to  transport  more 
sediment  (due  to  the  reduction  in  wave  height  and  angle  in  the  diffraction 
shadow  zone).  The  remaining  difference  is  also  due  to  the  volume  of  sediment 
being  impounded;  i.e.,  the  distance  and  extent  of  change  the  presence  of  the 
groin  causes  upcoast  and  downcoast. 

b.  Comparison  of  Cases  4.2c  and  4. 2d.  The  variations  between  cases  4.2c  and 
4. 2d  are  very  similar  to  the  differences  between  cases  4.2a  and  4.2b  as  would 

be  expected  with  a  groin  field  (here,  three  groins)  as  compared  with  a  single 
groin  (see  Figs.  9  and  10).  There  is,  however,  one  additional  feature  which 
can  be  attributed  to  the  additional  groins.  Note  that  in  the  direction  of 
littoral  drift,  the  size  of  the  fillet  Is  decreasing.  This  is  due  to  the 
updrift  beach  having  an  uninterrupted  supply  of  sediment  while  the  downdrift 
groin  compartments  are  supplied  sand  at  a  rate  determined  by  the  bypassing. 

Part  of  this  feature  may  also  be  due  to  the  system  not  having  attained 
complete  equilibrium. 

The  effects  of  the  fixed  boundary  conditions  are  evident  on  all  cases 
run.  In  these  example  cases,  the  boundaries  are  clearly  too  close  to  the 
structure  to  provide  a  proper  representation  of  the  fillet  contours. 

3 .  Simulations  of  Sediment  Transport  of  Dredge  Disposal  in  the  Vicinity  of 

Oregon  Inlet. 

Hypothetical  dredge  disposal  movement  in  the  nearshore  but  beyond  what 
is  normally  the  surf  zone  at  Oregon  Inlet's  adjacent  beach  to  the  south  was 
modeled.  In  order  to  do  these  simulations,  the  program  was  altered  such  that 
for  every  n 1?)  iteration  (time  periods),  the  contours  were  shifted  seaward 
to  simulate  the  addition  of  dredged  sediment  disposal.  The  program  presented 
in  Appendix  B  does  require  slight  modification  to  simulate  this  situation. 

In  general,  the  fifth  and  sixth  contours  were  shifted  seaward  on  a 
monthly  basis  to  simulate  the  disposal  of  121,000  cubic  yards  of  sediment. 
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8  contours  not  shown 


Longshore  Distance,  < 

Fiqure  9.  Eouilibrium  planform,  case  4.2c 


In  all  these  simulations,  the  following  variables  were  held  constant:  (a)  a 
time- step  of  3  hours,  { b )  a  shoreline  length  of  10,000  feet,  (c)  a  longshore 
space-step  of  200  feet,  (d)  an  A  value  of  0.15  foot1/^  for  the  equilibrium 
profile  (see  Fia.  11),  (e)  a  berm  height  of  5.3  feet  with  a  beach  face  slope 
of  0.05,  and  (f)  a  duration  of  1  year.  The  wave  climate  was  provided  by 
the  U.S.  Army  Engineer  Waterways  Experiment  Station  Wave  Information  Study 
(WIS)  1975  data  and  was  initiated  at  different  times  of  the  year  as  indi¬ 
cated  in  the  specific  cases  below.  All  simulations,  prior  to  any  addition 
of  sediment,  used  the  bathymetry  shown  in  Figure  12.  The  shoreline 
(relative  to  mean  low  water,  MLW)  was  scaled  from  a  bathymetry-topography 
survey  provided  by  the  U.  S.  Army  Engineer  District,  Wilmington.  The  initial 
offshore  bathymetry  was  computed  according  to  the  equilibrium  profile  and  the 
0-foot  contour;  i.e.,  the  profile  was  shifted  seaward  or  landward, 
accordingly,  (see  App.  C.)  The  boundary  profiles  were  fixed  throughout  the 
simulations.  The  variation  of  C0FF  outside  the  surf  zone  was  used  because  of 
the  importance  of  the  time  rate  of  change  in  this  simulation.  Table  1 
presents  the  percentage  of  sediment  which  moves  out  of  the  control  volume 
(i.e.,  imaginary  boundaries  around  the  area  where  sediment  was  added) 
directly  onshore  and  the  percentage  of  sediment  remaining  in  the  control 
volume  at  the  conclusion  of  the  simulation  for  each  of  the  cases.  In 
addition,  a  seventh  (case  3)  and  eighth  (case  4)  were  modeled.  In  Case  3, 
the  only  difference  was  that  sediment  was  placed  at  the  11-  and  14-foot 
contours.  Case  4,  however,  was  quite  different  and  will  be  described  in 
detail  later.  It  has  a  20,000-foot  shoreline,  a  longshore  space-step  of  400 
feet,  and  sediment  was  added  on  a  weekly  basis.  Also,  the  resolution  in  the 
profile  was  better. 

a .  Specific  Cases. 

(1)  Case  2. a.  In  order  to  provide  insight  for  the  interpretation  of  the 
other  modeling  efforts,  a  simulation  of  the  shoreline  evolution  using  the 
January  to  December  WIS  time  series,  with  no  addition  of  sediment,  was 
carried  out.  As  expected,  the  contours  almost  attain  an  equilibrium  planform 
shape  (i.e.,  straight  and  parallel  between  the  fixed  end  profiles;  they  do 
not,  however,  become  aligned  parallel  to  the  base  line  because  of  the  end 
conditions).  Because  of  the  scales  involved,  alongshore  versus 
onshore-offshore,  plotting  the  contours  without  distortion  does  not  yield 
much  information.  Appendix  C  provides  a  listing  of  the  final  contours  for 
all  the  cases  modeled. 

(2)  Case  2.b.  The  only  difference  between  cases  2. a  and  2.b  is  the 
suppression  of  the  WIS  wave  angle  which  was  set  equal  to  zero  (i.e.,  wave 
crest  approach  is  shore-parallel  at  the  offshore  boundary  of  the  model). 

This  does  not  cause  the  longshore  sediment  transport  to  vanish  completely. 
There  are  still  local  gradients  in  the  contours  which  cause  refraction  and 
relative  angles  between  wave  crest  and  contour,  thereby  driving  the  longshore 
sediment  transport  (even  if  refraction  was  not  considered,  the  local  angle 
between  the  wave  crest  and  contour  would  cause  sediment  transport).  Note  the 
larger  onshore  transport  (Table  1)  for  this  case  compared  with  Case  2. a. 

This  is  due  to  the  reduction  in  longshore  transport  caused  by  the  wave  angle 
of  0°.  The  model  still  tries  to  smooth  the  contour  lines;  however,  more  of 
the  smoothing  for  the  present  case  must  be  done  by  onshore-offshore  transport. 


Stepped  version  of  equilibrium  profile  used  in  the  Oregon  Inlet  model inq. 


imulations.  (The  scale  for  case  4  was  twice  the  scale  shown. 
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Table  1.  Summary  of  results  at  Oregon  Inlet. 


Case 

No. 

Description 

Pet  Onshore  out  of 
control  volume 

Pet  Remaining 
in  control  volume 

2. a 

No  sediment  added, 

WIS  waves 

Jan.  -  Dec. 

Onshore  Movement 
(992  yd3) 

Increase 
(14,148  yd3) 

2.b 

No  sediment  added. 

WIS  waves  (a  =0°' 

Jan.  -  Dec. 

Onshore  Movement 
(1624  yd3) 

Increase 
(9,356  yd3) 

2. cl 

121,000  yd3  added 
monthly,  WIS  waves 

Jan  -  Dec. 

31 .7 

(460,264  yd3) 

38.6 

(559,984  yd3) 

2.c2 

121,000  yd3  added 
monthly,  WIS  waves 
Apr.  -  Mar. 

32.1 

(466,160  yd3) 

36.9 

(535,392  yd3) 

2.c3 

121,000  yd3  added 
monthly,  WIS  waves 
July  -  June. 

28.6 

(415,784  yd3) 

47.0 

(682,088  yd3) 

** 

O 

_ 

121,000  yd3  added 
monthly,  WIS  waves 
Oct.  -  Sept. 

27.2 

(395,556  yd3) 

46.8 

(670,848  yd3) 

3 

121,000  yd3  added 
monthly  at  the  11- 
and  14-foot  contours 
WIS  waves,  Jan.  -  Dec. 

8.9  * 

(32,164  yd3) 

78.0 

(283,016  yd3) 

4 

27,923  yd3  added 
weekly  on  the  7- 
8-,  9-,  and  10-foot 
contours,  WIS  waves 
Jan.  -  Dec. 

19.0 

(275,796  yd3) 

47.4 

(687,525  yd3) 

*  After  17  weeks,  the  addition  of  sand  caused  contours  to  cross.  Prior 
sediment  added  was  363,000  yd3.  Problem  was  rectified;  however,  case  was 
not  rerun. 
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(3)  Case  2. cl.  In  this  simulation,  sediment  is  added  to  the  system  each 
month.  It  was  simulated  by  advancing  the  7-  and  11-foot  contours  on  a 
monthly  basis  to  represent  121,000  cubic  yards  per  month.  Specifically,  the 
sand  volumes  were  "tapered"  starting  at  the  center  of  the  nourished  area  over 
a  distance  of  +  2,700  feet  from  the  center.  Table  2  presents  the  monthly  by 
values  for  the  blocks  between  the  7-  to  11-foot  contours  and  the  11-  to-14 
foot  contours.  Figure  13  shows  the  planform  ty  values  added  monthly.  WIS 
waves  were  used  with  the  sequence  being  the  normal  calendar  year,  January 
through  December. 


contour 

depth 

increment 


width  of 
grid 


OFFSHORE 

DIRECTION 


Figure  13.  Monthly  incremental  values  of  _\y  due  to  dredge  disposal 

illustrated  for  the  block  between  7-  and  11-foe1;  contours. 


The  initial  and  final  fifth  and  sixth  contours  have  been  plotted  in 
Figures  14  and  15.  The  first  figure  has  no  distortion;  the  second  is 
distorted  10  to  1.  The  simulation  predicts  that  31.7  percent  of  the  dredge 
disposal  will  move  shoreward  out  of  the  control  volume.  An  additional  29.7 
percent  efflux  occurs  in  the  offshore  and  longshore  directions,  leaving  only 
38.6  percent  of  the  total  amount  of  sediment  added  remaining  in  the  control 
volume.  It  is  not  clear  what  quantity  of  the  sediment  leaving  in  the 
longshore  direction  would  reach  shore.  It  is  conceivable  that  most  of  this 
sediment  would  eventually  reach  the  surf  zone.  The  rate  at  which  this 
material  would  move  ashore  would  be  expected  to  be  slower  than  the  rate  at 
which  the  large  mounds  would  move  ashore  because  the  deviation  of  the  profile 
from  equilibrium  is  much  less. 

(4)  Cases  2.c2,  2.c3,  and  2.c4.  The  next  three  simulations  were  the 
same  as  2. cl  except  the  time  series  of  wave  events  has  been  seasonally 
altered.  Cases  2.c2,  2.c3  and  2.c4  use  the  1975  wave  climate  from  April 
through  March,  July  through  June,  and  October  through  September,  respectively. 
The  maximum  variation  is  about  5  percent  for  the  sediment  volume  moving 
onshore,  and  about  10  percent  for  the  volume  remaining.  The  variation  in  the 
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Table  2.  Monthly  values  of  Ay  for  the  steps  located  between  the  7-  to 
10-foot  contours  and  the  11-  to  14-foot  contours. 


Value  of  I 

Monthly  ay  valup  (ft) 

for  steps  between 

7-  and  11-foot  contours 

11-  and  14-foot  contours 

26 

145.8 

194.4 

25,27 

135.4 

180.5 

24,28 

125.0 

166.6 

23,29 

114.6 

152.7 

22,30 

104.1 

138.9 

21,31 

93.7 

125.0 

20,32 

83.3 

111.1 

19,33 

72.9 

97.2 

18,34 

62.5 

83.3 

17,35 

52.1 

69.4 

16,36 

41.7 

55.5 

15,37 

31.2 

41.7 

14,38 

20.8 

27.8 

13,39 

10.4 

13.9 

All  Others 

0 

0 
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Figure  14.  Initial  and  final  7-  and  11-foot  contours  (no  distortion). 
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al  and  final  contours  for  case  2. cl  [y(I,5)  and  y(I,6)] 


quantity  moving  onshore  could  be  caused  by  waves  that  first  tend  to  move  more 
sediment  longshore;  then,  the  waves  that  transport  more  sediment  onshore  have 
a  less  out-of-equil ibrium  profile  to  cause  movement  upon.  The  variation  in 
percentage  remaining  is  due  to  the  variation  of  the  time  series  of  the  wave 
climate,  with  the  last  month  in  the  simulation  being  especially  important. 

(5)  Case  3.  Instead  of  extending  the  7-  and  11-foot  contours  monthly  to 
simulate  the  disposal  of  dredged  sediments,  the  11-  and  14-foot  contours  were 
extended  (194.4  feet  each  at  the  center  of  the  disposal  area).  This  case  was 
modeled  because  the  larger  available  dredge  could  not  dump  in  more  shallow 
water.  The  reduction  and  increase  in  the  percent  of  onshore  volume  and  the 
percent  volume  remaining  (8.9  percent  and  78.0  percent,  respectively) 
demonstrate  the  sensitivity  of  the  depths  investigated.  Qualitatively,  these 
depths  are  the  depths  to  which  offshore  bars  occur  along  the  Atlantic  U.S. 
coast. 

(6)  Case  4.  Further  investigation  of  the  disposal  process  demonstrated 
the  need  for  an  11,000-foot  shore-parallel  disposal  length  with  the  sediment 
placed  at  the  11-foot  contour  building  to  about  7  feet.  It  was  decided  to 
model  this  physical  situation  also.  The  total  shoreline  length  was  changed 
to  20,000  feet,  and  the  space  step  to  400  feet;  the  length  of  the  disposal 
area  in  the  longshore  direction  was  increased  to  10,800  feet.  The  resolution 
in  the  vicinity  of  the  depths  of  the  dump  was  improved  by  adding  the 
additional  contours  and  the  profile  is  shown  in  Figure  16.  As  in  the  other 
seven  cases,  1,452,000  cubic  yards  was  added  annually  to  the  system;  however, 
the  addition  was  accomplished  on  a  weekly  basis  (27,923  cubic  yards  per 
week).  Sediment  was  still  added  by  extending  the  contours  seaward,  but 
rather  than  placing  one-fourth  of  the  sediment  at  each  of  the  four  contours, 
the  volumes  were  determined  based  on  the  trapezoidal  cross  section  shown  in 
Figure  17.  This  cross  section  more  closely  resembles  the  disposal  mound 
formed  by  hopper  dredging.  The  incremental  values  Figure  18  show,  in 
planform,  the  extension  of  the  contours  to  simulate  the  weekly  sediment 
addition  at  the  8-foot  contour. 

A  schematic  illustration  of  the  sediment  transported  from  the  nourished 
region  is  presented  in  Appendix  C.  Nineteen  percent  of  the  sediment  added 
moved  directly  onshore  out  of  the  control  volume. 

b .  Conclusions  for  the  Movement  of  Disposed  Sediment  i n  the  Vi ci ni ty  of 
Oregon  Irnei.  The  computer  simulations,  tempered  with  engineering  judgment, 

demonstrate  that  between  15  and  35  percent  of  the  material  added  to  the  7- 
and  11-foot  contours, or  to  the  7-  8-  9-,  and  10-foot  contours  would  be 
transported  into  the  nearshore  transport  system  during  the  first  year.  If 
the  disposal  process  was  continued,  the  system  would  approach  steady  state  in 
terms  of  the  volume  of  deposited  material  residing  offshore. 

For  the  case  of  sediment  addition  at  the  11-  and  14-foot  contours,  the 
computer  simulations,  tempered  with  engineering  judgment,  show  that  between  5 
and  25  percent  of  the  material  added  would  be  transported  into  the  nearshore 
transport  system  during  the  first  year. 
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Figure  16.  Stepped  version  of  equilibrium  profile  used  in  the  Oregon  Inlet  modeling,  h=Ay 
(A=0.15  feet1/3),  case  4.  Note  the  resolution  at  7,  8,  9,  and  10. feet. 


Figure  17.  Shore-perpendicular  cross  section  of  disposal 
mound.  The  volumes  represent  the  volume  per¬ 
centage  of  the  trapezoidal  section  between  contours 
and  therefore,  the  quantity  of  sediment  added  to 
the  7-,  8-,  9-,  and  10-  foot  contours. 


Figure  18.  Incremental  values  of  Ay  due  to  dredge  disposal, 
illustrated  for  the  block  between  8-  and  9-foot 
contours  (case  4). 


4 .  Simulation  of  the  Longs ho r e_ Sand  Transport  Study  at  Channel  Islands  Harbor , 
Cal i fornia . 

The  CIH  Longshore  Sand  Transport  Study  (Bruno,  et  al . ,  1981)  was  the 
only  field  study  found  suitable  for  verification  purposes.  Wave  data 
collected  included  the  LEO  data  and  a  two  pressure-sensor  gage  array. 

Although  the  pressure  gages  were  not  in  operation  throughout  the  study,  it 
was  expected  that  the  data  they  produced  would  be  superior  to  that  of  the  LEO 
data.  However,  these  data  were  not  available  in  a  reduced  form,  so  the  LEO 
data  were  used.  An  adjustment  of  11°  was  made  to  the  breaker  angle  to  orient 
the  angle  with  respect  to  the  base  line,  rather  than  to  the  local  shoreline 
orientation  angle.  Observations  had  been  taken  twice  daily  at  three 
locations;  the  middle  location  was  used  (observer  No.  5714).  Waves  which 
approached  the  shoreline  at  angles  too  large  to  have  originated  in  a  depth  of 
10  meters,  according  to  Snell's  law,  were  set  equal  to  90°  at  that  depth 
(crest  of  wave  perpendicular  to  the  baseline).  The  10-meter  depth  was  chosen 
because  it  is  the  approximate  depth  at  the  tip  of  the  offshore  breakwater 
(for  this  reason,  it  was  also  chosen  as  the  depth  of  the  step  beyond  the  y(I, 
JMAX  +  2)th  contour).  It  was  assumed  that  each  of  the  two  daily 
observations  occurred  for  12  hours  and  using  a  time-step  of  6  hours,  this 
meant  two  time-steps  per  wave.  In  cases  where  parts  of  the  wave  data  (Hb, 
at,,  or  T)  were  missed  by  the  observer  or  were  equal  to  zero,  the  data  were 
ignored  (no  computations  were  made),  but  the  time  was  included.  Because  the 
time  rate  of  change  is  important  for  this  simulation,  the  variation  of  Coff 
outside  the  break  point  was  used. 

The  period  chosen  to  model  was  20  April  through  1  December  1976.  The 
initial  survey  was  taken  after  dredging  of  the  sediment  trap  and  for  this 
reason  was  known  to  be  out.  of  equilibrium.  The  bathymetric  surveys  were  con¬ 
ducted  using  several  methods,  the  most  advanced  being  a  Lighter  Amphibious 
Resupply  Cargo  vessel  (LARC)  proceeding  along  shore-perpendicular  lines 
(approximately  in  the  vicinity  of  each  survey  station)  taking  fathometer 
readings  every  10  seconds,  with  positioning  systems  trilaterating  the 
vessel's  position  concurrently.  These  data  were  recorded  on  tape.  The 
beach-face  data  were  taken  using  standard  surveying  methods.  Because  the 
data  fluctuated  randomly  about  the  stations,  depending  on  the  speed  of 
the  craft,  the  (x,  y)  coordinate  positions  had  to  be  altered  to  fixed  changes 
in  x  and  y.  This  was  accomplished  using  an  interpolation  routine.  The  x 
values  were  made  to  coincide  with  the  stations  used  in  the  surveys,  and  the  y 
values  were  determined  at  100-foot  intervals  beginning  from  the  base  line. 
Stations  100+00  and  118+00  were  located  at  the  north  jetty  and  termination  of 
the  detached  breakwater,  respectively  (these  correspond  to  I  values  of  16.5 
and  34.5  in  the  model).  See  Figure  19. 

Monotonic  profiles  of  the  form  h  =  A(y  -  ydel )2/3  were  fit  to  the  data 
along  each  station  line,  "ydel"  represents  the  zero  location  of  the  fitted 
shoreline,  the  value  of  which  was  unknown.  Because  dredging  had  been  done  in 
the  lee  of  the  breakwater,  there  was  no  reason  to  expect  the  A  value  to 
correspond  to  the  value  upcoast  where  the  influence  of  the  structure  and  the 
dredging  was  negligible.  For  this  reason,  the  profiles  of  Stations  122+00 
through  134+00  were  evaluated  separately  to  determine  an  A  value  for  the 
equilibrium  profile  to  be  used  in  the  numerical  model.  For  the  detailed 
method  used  (LaGrange  Multipliers  and  Newton-Raphson  Method  for  nonlinear 
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Figure  19.  Idealized  numerical  model  representation  of  offshore 
breakwater  at  Channel  Islands  Harbor,  California. 
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equations)  and  the  computer  programs  see  Appendix  D.  The  two  values  obtained 
for  the  surveys  of  20  April  and  1  December  1976  were  averaged  to  obtain  the 
value  used  in  the  model,  A  =  0.2606.  Stations  101+00  through  121+00  were 
treated  separately  for  the  purpose  of  obtaining  values  with  which  to  initial¬ 
ize  those  parts  of  the  contours  in  the  model  and  for  comparison  of  the  model 
predictions  with  the  prototype  values.  Note  that  although  the  breakwater 
extends  only  to  about  Station  118+00,  the  influence  of  the  structure  and 
dredging  extends  beyond  that  location  and  so,  although  arbitrary,  the  121+00 
station  was  chosen  as  the  dividing  line.  The  initial  and  final  values  of  the 
scaling  parameter  A  for  the  profiles  were  0.3233  and  0.3528,  respectively. 
Because  the  initial  shoreline  is  so  irregular,  a  discontinuity  between  121+00 
and  122+00  is  not  evident. 

One  further  idealization  was  made.  The  jetty-breakwater  system  was 
idealized  as  shown  in  Figure  19.  This  was  required  to  simplify  the  physical 
situation,  and  although  waves,  currents,  and  sediment  do  pass  through  the 
opening  in  the  prototype,  it  is  hoped  that  they  are  of  secondary  importance. 

The  results  of  the  numerical  modeling  of  Channel  Islands  Harbor  are 
presented  in  Figures  20  and  21.  The  first  figure  presents  the  shoreline 
contour  (depth  =  0);  the  second  figure  presents  the  farthest  offshore, 
modeled  contour.  In  both  cases,  the  initial  shoreline  represents  the  model 
and  prototype  (after  fitting  of  the  profiles).  The  initial  shoreline  contour 
is  further  offshore  along  the  section  of  beach  beyond  the  end  of  the 
breakwater,  while  in  the  lee  of  the  breakwater,  as  would  be  expected  after 
dredging,  the  shoreline  is  closer  to  the  base  line.  The  final  prototype 
contour  has  undergone  erosion  along  the  reach  beyond  the  tip  of  the 
structure,  and  accretion  in  the  lee. 

The  model's  shoreline  contour  has  undergone  similar  changes,  and  on  the 
average,  represents  the  final  prototype  contour  quite  well.  The  JMAXl!) 
contour  has  been  displaced  quite  similarly  to  the  shoreline  contour  with 
shoreward  movement  (erosion)  along  the  reach  beyond  the  tip  of  the  breakwater 
and  seaward  movement  (accretion)  within.  It  appears  that  the  final  model’s 
shoreline  has  predicted  too  much  erosion  and  not  enough  accretion.  Several 
parameters  could  be  incorrect,  with  the  onshore-offshore  sediment  transport 
rate  coefficient,  Cqpp,  perhaps  the  most  likely.  Overall,  the  model  seemed 
to  predict  reasonable  values  oi  the  contours. 


V.  SUMMARY  AND  RECOMMENDATIONS 


Some  of  the  parameters  that  the  model  not  include  are  important  and 
should  be  mentioned.  As  stated  previa  Jt  the  model  does  not  include  bar 
formation.  This  is  precluded  by  an  n-line  system.  There  are  no  provisions 
for  water  level  fluctuations  or  currents.  Improvement  to  the  model  could 
also  be  facilitated  with  better  longshore  and  cross-shore  sediment  transport 
relationships.  A  more  reliable  equation  for  distribution  of  sediment 
transport  across  the  surf  zone  would  also  be  helpful  (or  further  testing  and 
calibration  of  the  equation  proposed  herein).  Finally,  combining  refraction 
and  diffraction  using  equations  to  predict  their  combined  effect  would 
improve  the  wave  field.  The  program  was  constructed  such  that  improvement 
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Figure  21.  CIH  simulation  of  (JMAX)^  contour,  20  April 
1  December  1976  (from  LEO  data). 
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could  be  accomplished  with  minimum  effort.  Therefore,  if  a  more  suitable 
equation  becomes  available,  the  change  of  a  subroutine  should  be  sufficient 
for  implementation  of  the  equation. 

Although  the  model  is  limited  by  the  omission  of  the  aforementioned 
parameters,  it  is  reasonably  correct.  The  ability  to  simulate  various 
physical  situations  (shore-perpendicular  structures,  beach  fills,  breakwater 
and  shore-perpendicular  structures)  has  been  demonstrated.  In  the  CIH 
simulation  where  the  data  were  first  transformed  to  monotonically  decreasing 
contours  and  LEO  wave  data  were  used,  the  model  still  predicts  the  prototype 
shoreline  changes  in  a  reasonable  fashion. 

Further  research  and  model  development  should  include  exercising  the 
model  in  a  number  of  different  situations.  Several  theoretical  cases  should 
be  simulated,  which  if  analyzed  properly,  would  provide  a  tool  for  the 
coastal  engineer.  Combined  refraction  and  diffraction  should  be  included,  if 
possible,  along  with  any  of  the  aforementioned  parameters  which  have  been 
omitted  and  for  which  relationships  exist.  Perhaps  the  most  difficult  prob¬ 
lem  to  researchers  working  on  modeling  sediment  transport  in  the  vicinity  of 
structures  is  the  availability  of  field  data.  High-quality  concurrent  wave 
and  bathymetric  change  data  in  the  vicinity  of  coastal  structures  do  not 
exist.  One  suggested  field  experiment  is  to  monitor  changes  both  updrift  and 
downdrift  of  a  jettied  inlet  which  has  a  bypassing  plant.  Monitoring  should 
begin  immediately  after  bypassing,  when  the  profiles  are  out  of  equilibrium. 
The  recorded  bathymetric  and  wave  data  would  then  provide  data  with  which  to 
calibrate,  verify,  and  evaluate  the  existing  models. 
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APPENDIX  A 


DISCUSSION  OF  CONSTANTS  AND  SOME  OF  THE  VARIABLES 
REQUIRED  BY  THE  MODEL 


Establishing  the  grid-contour  system  requires  several  variables.  IMAX 
represents  the  number  of  cross-shore  grid  lines  desired  and  JMAX  the  number 
of  contours  simulated.  DX  represents  the  spacing  between  the  IMAX  grid  lines 
and  DY  the  spacing  between  the  contours.  DX  is  a  value  which  must  be  chosen 
along  with  IMAX  and  JMAX  such  that  sufficient  detail  is  obtained  where 
necessary  (e.g.,  in  the  shadow  zone,  if  diffraction  effects  are  believed  to 
be  very  important,  DX  must  be  assigned  a  sufficiently  small  value  so  that  at 
least  some  points  lie  within  the  shadow  zone  for  the  larger  wave  angles).  DY 
is  not  a  constant,  but  a  dimensional  array  which  is  computed  by  the  model 
according  to  the  contour  location.  Once  the  depths  of  contours  to  be  modeled 
are  chosen,  the  initialization  of  DY  and  the  y  values  are  computed  with  the 
following  equation  a^ter  Dean,  1977 

h  =  A  y2/3  (A-l ) 


where  h  is  the  depth,  y  is  the  offshore  distance  and  A  is  the  scaling  parameter 
Dean  gives  values  for  A  for  several  diameter  sediments;  however,  if  long-term 
beach  profiles  are  available  for  the  site  being  modeled,  the  modeler  may  want 
to  choose  a  slightly  different  A  value  to  more  closely  match  the  site-specific 
beach  profile.  Figure  A-l  presents  values  of  A  versus  diameter  (after  Moore, 
1982).  The  model  is  programmed  to  input  the  h(I,J)  values  (depths  as  shown  in 
Figure  1,  called  DEEP  (I,J)  in  the  program)  read  in  the  value  of  A  (called 
ADEAN  in  program)  and  it  then  computes  the  y  values.  Also  shown  in  Figure  1 
is  the  height  of  the  berm  (BERM)  and  this  value,  along  with  the  beach-face 
slope  (SFACE),  is  required  as  program  input  and  can  be  obtained  from  beach 
profile  site  data.  Because  the  model  does  not  include  water  level 
fluctuations  such  as  tides,  all  values  are  to  be  referenced  to  a  chosen 
datum.  Other  geometrical  constants  depending  on  the  site  include  SJETTY  (the 
length  of  the  jetty),  MMAX  (the  number  of  structures  to  be  input),  and  IJET 
(M),  M  =  1,2, ...MMAX  (the  smaller  I  value  adjacent  to  the  Mth  structure's 
location).  If  no  structure  is  required,  as  in  a  beach  fill,  the  value  of 
SJETTY  must  be  entered  as  0.0,  with  MMAX  and  IJET  (M)  entered  as  1  and 
(IMAX/2),  respectively.  As  set  up  presently,  the  groin  locations  must  be 
equally  spaced. 

One  constant  used  throughout  the  program  is  the  breaking  wave  criteria 
(CAPPA  in  the  program)  equal  to  0.78.  It  is  required  in  several  different 
computations  and  always  governs  the  maximum  wave  height  allowed  according  to 
the  depth. 

Another  group  of  variables  assigned  values  within  the  program  is  the 
sediment  and  fluid  properties.  These  include  fluid  mass  density,  sediment 
mass  density,  porosity,  and  the  angle  of  repose  (e.g.,  RHO  =  1.99,  RHOS  = 

5.14,  POROS  =  0.40,  and  REPOSE  =  32°,  respectively).  The  values  can  easily  be 
changed  to  reflect  site  conditions. 
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Figure  A-l.  A  versus  sediment  diameter  (after  Moore,  1982) 


Another  very  important  set  of  constants  is  the  constant  chosen  for  the 
longshore  and  cross-shore  components  of  sediment  transport.  Equation  (27), 
the  total  longshore  transport  equation,  contains  the  constant  C'  equal  to 


where  K  =  0.77  (Komar  and  Inman,  1970) 

g  is  the  acceleration  of  gravity  (32.17  ft/sec2) 

o s  and  n  are  the  mass  densities  of  the  sediment 
and  the  seawater  (5.14  and  1.99  slugs  per  cubic  feet, 
respectively 

p  is  the  porosity  (0.40),  and 
<  is  taken  as  0.78. 

Using  these  values  to  compute  C‘  (TKSI  in  the  program),  a  value  of  0.325  is 
obtained.  It  is  stressed  that  if  any  of  these  values  are  different  for  the 
site  to  be  modeled,  they  should  be  changed  and  the  program  will  compute 
another  value  for  C1. 

The  parameter  CoFF  is  an  "activity  factor"  which,  based  on  earlier  work 
primarily  within  the  surf  zone,  was  found  to  be 

CQFF  =  10'5  ft/s,  h  <  hb 

To  generalize  this  concept  for  transport  seaward  of  the  surf  zone,  the 
wave  energy  dissipation  per  unit  volume  was  utilized  as  a  measure  of 
mobilization  of  the  bottom  sediment.  Inside  the  surf  zone,  the  dominant  wave 
energy  dissipation  is  caused  by  wave  breaking;  outside  the  surf  zone,  the 
dominant  mode  of  wave  energy  dissipation  is  due  to  bottom  friction.  These 
two  components  will  be  denoted  by  and  D?,  respectively. 

(a)  Energy  Dissipation  by  Wave  Breaking.  The  wave  energy  dissipation 
per  unit  volume  by  wave  breaking,  Dj,  is 


(E  CQ) 


(A-3) 


which,  employing  the  spilling  breaker  assumption  (H  =  <h)  within  the  surf 
zone,  can  be  shown  to  be 
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(A-4) 


n  _  5  3/2  2.1/2  ah 

°1  =  IS  pg  K  h  55 


5  3/2k2a3/2 

25  ^9  *  A 


(A  5) 


in  which  A  is  the  scale  parameter  in  the  equilibrium  beach  profile 


h  (y)  =  Ay'1 


(A-6) 


( b )  Energy  Dissipation  by  Bottom  Friction.  The  wave  energy  dissipation 
per  unit  volume  due  to  bottom  friction,  D,,,  is 


°2  =  h  T  ub  =  fi  pCf  lutJ  ub 


(A-7 ) 


in  which  Cf  is  a  bottom  friction  coefficient,  is  the  bottom  water 
particle  velocity  and  the  overbar  indicates  a  time  average.  For  linear 
waves,  equation  (A-7)  can  be  reduced  to 


n  _  1  p  r  HV 


2  ~  6ir  fi  f  "7 


(A-8) 


sinh  kh 


The  activity  coefficient  CnFF,  outside  the  surf  zone,  is  expressed  as 


'OFF  *  f  Sf  x  10'5  ft/5’  h  >  hb 


(A-9) 


'OFF  =  5T 


% 

'fa  I 


H  1 
sinh  kh, 


x  10' 


( A- 10 ) 


in  which  r  is  a  parameter  relating  the  efficiency  with  which  breaking  wave 
energy  (which  occurs  primarily  near  the  water  surface)  mobilizes  the  sediment 
bottom  (0  <  r  _<  1).  Herein,  r  is  taken  to  be  one. 

Figure  A-2  presents  an  example  of  the  variation  of  the  activity 
coefficient  versus  relative  depth  for  a  particular  wave  period  and  deep  water 
wave  height.  It  is  seen  that  the  activity  coefficient  reduces  rapidly  with 
increasing  depth. 

The  value  of  Coff  f°r  the  physical  modeling  of  Savage's  (1959)  data 
was  taken  as  10"*  feet  per  second.  Perl  in  (1978)  presents  some  rationale 
for  choosing  a  value  of  Coff;  however,  very  little  testing  has  been  done 
and  none  is  based  on  actual  field  measurement. 
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Finally,  wave  data  are  read  into  the  program  and  the  simulation  begins. 
(For  information  regarding  "Read  Formats"  for  the  various  constants  and 
variables,  see  Appendix  E).  Wave  data  required  are  wave  height,  wave  period, 
wave  angle  relative  to  the  x-axis  of  the  model  at  a  depth,  WDEPTH  and  the 
duration  of  the  wave  climate  (HS,  T,  ALPWIS,  and  a  combination  of  NTIMES  x 
DELT,  respectively,  in  the  model).  As  is  always  the  case  with  numerical 
models,  the  time  step  and  space  steps  are  very  important  to  both  stability 
and  accuracy.  Time  steps  on  the  order  of  3  to  6  hours  (10,800  to  21,600 
seconds)  or  less  are  recommended.  However,  the  complexity  of  the  bathymetry, 
variation  and  time  series  of  the  wave  data,  constants  used  (especially 
Coff)  along  with  several  other  factors,  greatly  influences  the  stability 
and  accuracy  of  the  results. 

Table  A-l  lists  several  of  the  important  variables  in  the  computer 
program. 


Table  A-l.  List  of  important  variables  in  the  program 


ABAND  The  input  banded  matrix  which  stores  the  values  from  equation  (37) 

ADEAN  The  value  of  the  scaling  parameter  in  the  equilibrium  beach  profile 

ALPHAS  The  angle  a  contour  makes  with  the  x-direction  base  line 

(counter-clockwise  is  positive) 

ALPWIS  The  angle  (-90°  to  *-90°)  the  wave  crest  makes  with  the 
x-direction  (counter-clockwis'-  is  positive) 

AMP  The  amplitude  of  the  diffracted  wave  in  the  shadow  zone 

ANGGEN  The  wave  angle  at  a  depth,  WDEPTH 

ANGLOC  The  local  contour  orientation  angle 

AWARE  See  equations  (36)  and  (37) 

BERM  The  height  of  the  berm  above  water  level 

BMATRX  The  matrix  which,  upon  solution  of  the  banded  matrix  problem 

yields  the  new  y  values 

C  The  wave  celerity 

CAPPA  The  breaking  wave  index 

CC  Constant  which  establishes  the  width  of  the  distribution  of 

sediment  transport  across  the  surf  zone 

CG  The  group  velocity  throughout  the  wave  field 

CGEN  The  linear  wave  theory  celerity  at  a  depth,  WDEPTH 
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CGGEN 

CO 

COFF 

CONST 

C0NST6 

DEEP 

DEEPB 

DEEPBI 

DELT 

DIAM 

DISTR 

DX 

DY 

ELO 

ELTIP 

EPS 

G 

GAMMA 

H 

HB 

HBI 

HBQ 


The  linear  wave  theory  group  velocity  at  a  depth,  WDEPTH 

The  deepwater,  linear  wave  theory  wave  celerity 

The  onshore-offshore  transport  rate  coefficient  within  the  surf 
zone 

The  constant  in  the  longshore  sediment  transport  relationship 
(0.77) 

The  space  step,  DX,  multiplied  by  the  activity  coefficient 
The  water  depth  at  any  grid  location 

The  initial  breaking  depth  along  each  profile  (between  adjacent 
profiles) 

The  initial  breaking  depth  along  each  profile  (at  each  profile, 
rather  than  between  them) 

The  time-step  in  seconds  (DELT  x  NTIMES  =  wave  condition  duration) 
The  mean  diameter  of  the  sediment  particles 
See  equations  (36)  and  (37) 

The  alongshore  space-step  in  the  x-direction  (distance  between  I 
values) 

The  onshore-offshore  space-step  in  the  y-direction  as  defined  by 
the  stepped  profile 

The  deepwater,  linear  wave  theory  wavelength 

The  wavelength  at  the  tip  of  the  structure 

The  change  in  the  wave  number  which  is  acceptably  small 

The  acceleration  of  gravity  (32.17  feet/second^) 

The  specific  weight  of  seawater 

The  wave  height  throughout  the  wave  field 

The  maximum  wave  height  which  could  exist  throughout  the  wave 
field  (where  H  =  0.78  *  h) 

The  initial  breaking  wave  height  along  any  profile  at  the  y  values 
rather  than  between  them 

The  initial  breaking  wave  height  along  any  profile,  between 
adjacent  profiles 
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HGEN 


HS 

I 

I  BREAK 
IJET 

I  MAX 
J 

JMAX 

JUSE 

J1 

J2 

J1REF 

J2REF 

MMAX 

NITER 
NT  I  ME 
NTIMES 

NUN  IV 
PI 

POROS 

QX 

QXTOT 

QY 

R 


Average  wave  height  at  a  depth,  WDEPTH 
The  significant  wave  height  input 
The  longshore  grid  location 

The  leeward  side  of  the  initial  breaker  location  J  value 

Represents  the  lesser  I  value  adjacent  to  the  structure  (these 
must  be  evenly  spaced  alongshore) 

The  total  number  of  grid  points  in  the  x-direction  (alongshore) 

The  offshore  contour  location 

The  value  of  the  seawardmost  contour  simulated 

(JMAX  +  2)  the  seawardmost  contour  at  which  the  wave  field  is 
calculated 

Landward  contour  of  refraction  zone 

Seaward  contour  of  refraction  zone 

Landward  J  values  of  boundary  of  refraction  zone 

Seaward  J  values  of  boundary  of  refraction  zone 

The  number  of  shore-perpendicular  structures  to  be  simulated 
(present  maximum  of  16) 

The  counterindex  in  the  refraction  routine 

The  counterindex  in  the  time  simulation  "DO"  loop 

The  number  of  iterations  of  time-step,  DELT,  for  which  a 
particular  wave  Is  simulated 

The  total  number  of  time-steps  simulated  at  any  time 
The  value  of  it  =  3.141592654 
The  porosity  of  the  sediment 

The  longshore  sediment  transport  rate  at  a  specific  location 

The  total  alongshore  sediment  transport  rate  due  to  the  height  and 
angle  of  the  initial  breaking  wave 

The  onshore-offshore  sediment  transport  rate  at  a  specific  location 
See  equations  (36)  and  (37) 
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REPOSE 

RHO 

RHOND 

RHOS 

RK 

S3 

SFACE 
S JETTY 

SIGMA 

T 

TAU 

THETA 

THEATO 

TKSI 

TWOPI 

U 

UCRIT 

V 

WDEPTH 

WEQ 

XCOOR 

XDISTN 

Y 


The  angle  of  repose  of  the  sediment 
The  mass  density  of  seawater 

The  dimensionless  distance  from  the  tip  of  structure  where 
diffraction  is  initiated 

The  mass  density  of  sediment 

The  wave  number 

See  equations  (36)  and  (37) 

The  slope  of  the  shoreface 

The  length  of  the  shore-perpendicular  structure  (from  the  base 
line) 

The  wave  radian  frequency 
The  wave  period 

The  dissipative  interface  parameter 
The  wave  angle  throughout  the  wave  field 
The  wave  angle  at  the  tip  of  the  structure 
The  longshore  sediment  transport  rate  coefficient 
Twice  the  value  of  * 

See  equations  (36)  and  (37) 

The  critical  velocity  required  to  move  the  sediment  according  to 
the  She id's  diagram 

See  equations  (36)  and  (37) 

The  depth  of  water  in  meters  to  which  the  input  wave  conditions 
are  to  be  transformed 

The  equilibrium  profile  distance  between  contours  as  defined  by 
the  stepped  profile 

The  x-coordinate  where  the  wave  field  is  to  be  calculated. 

Together  with  YCOOR,  they  determine  whether  the  position  is  within 
or  beyond  the  diffraction  shadow  zone 

The  location  of  the  structure  along  the  shoreline  in  feet 

The  distance  offshore  to  the  contours 
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YCOOR 


1 

The  y-coordinate  where  the  wave  field  is  to  be  calculated. 

Together  with  XCOOR,  they  determine  whether  the  position  is  within 
or  beyond  the  diffraction  shadow  zone 

YDISS  The  value  of  y  after  the  use  of  the  dissipative  interface 

YOLD  The  previous  value  of  y 

YZERO  The  berm  contour  location 

Z1  See  equation  (37) 

12  See  equation  (37) 
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APPENDIX  B 
PROGRAM  LISTING 


i 


too  C*  . .  PROGRAM  IMPLICIT  SEOTRAN 

200  C-THIS  PROGRAM  IS  SET-UP  TO  HANDLE  MULTIPLE  GRO 1  NS  <  M' *  10  ) 

300  COMMON/ A/  C<60. 20) ,RK( GO. 20 ) . V ( 60. 201 .DEEP! 60. 20) . A l  PHAS1 GO. 20  I 

400  COMMON/ A A / Y  Z  E  RO ( 60 ) 

500  COMMON /BB/ WE 0( 60,20) 

600  COMMON/B/  THETA(60. 20) ,0XT0T(60( ,  0LDANGI6 0,20),  0060.201 

700  COMMON/C/  HI  60 . 20 ) . CGI  60 . 20 ) . HOLD ( 60 , 20 ) . HB ( 60 . 20 ) . VB ( 60 ) 

800  COMMON/N  USED/JUSE . T , CO . CGEN . CGGEN . ANGGEN . DX , BERM . THE  TAOI  101. MMA 

900  COMMON/D/ S I GMA.G.ELO, UMAX. I  MAX. PI .TW0P1  .PI02.HGEN.  I  JET  I  101 .SJETT 

lOOO  COMMON/F/AOEAN. REPOSE. DIAM 

1100  COMMON/ AAA/DELT .NTIMES 

1200  COMMON/COUNT/NUNI V 

1300  COMMON/EXPL/QVEXPIGO. 20) . V  IMP (60. 20) 

1400  DIMENSION  CHANGE  I  20 ) , HC I  10 ) . TC(  10 ) 

1500  DIMENSION  YOR I G I  60 . 20 )  . V ?E ROO I  60 ) . S ANGL E I  20 ) 

1600  NUNI V=0 

1700  JMAX=8 

1800  JUS  E  =  JMA  X  +■  2 

1900  I  MAX  =  50 

2000  PI =3  141592654 

2100  TWOP I -PI  * 2 . 

2200  P 1 02  =  P I /2  0 

2300  REPOSE  =  32 . *  TWOP 1/360 

2400  WRITEI6.732) 

2500  7  32  FORMAT!  . . . . . •  ) 

2600  WRITEI6.733) 

2700  733  FORMAT! 2X .' TO  WHAT  DEPTH  ARE  THE  WAVES  TO  BE  TRANSFORMED ' I 

2800  C'WDEPTH  MUST  BE  A  DEPTH  BEYOND  THE  END  OF  THE  STRUCT,  PREFERABLY  AT 
2900  C**DEEP( JMAX )  OR  GREATER! OR  ELSE  SNELL'S  LAW  DR  SHOAL  COULD  BLOWUP  IN 

3000  C**'DEEPER  WATER  ITS  IN  METERS  HERE' 

3100  READ (5,770)  WOE PTH 

3200  770  FORMAT! 10X , F 10  3 ) 

3300  WDEPTH=WDEPTH*3  28084 

3400  WR I TE ( 6 , 762  )  WOEPTH 

3500  762  F0RMAT(2X.  "THE  DEPTH  (IN  FT)  WAVES  TRANSFORMED  TO.  WDE PTH-  " 

3600  *  F10  3) 

3700  WR IT£(6,732) 

3800  WRITE(6,777) 

3900  777  FORMAT ( 2X I TS  TIME  FOR  SJETTY.  BERM.  ST  ACE .  AND  DIAM"./) 

4000  C*SJETTY  MUST  BE  MUCH  LESS  THAN  Y(l.JMAX) 

4100  READ! 5 . 776 )  S JE TT Y , BERM . SF ACE . D I  AM 

4200  776  FORMAT! 2F 10  3. F 10  4, F 10  3 ) 

4300  WR I TE ( 6 , 76 1 )  SJETTY 

4400  761  FORMAT! 2X THE  LENGTH  OF  THE  STRUCTURE.  SJETTY*  ' . F 10  3) 

4500  WR ITE(6.740)  BERM 

4600  740  FORMAT! 2X THE  HEIGHT  OF  THE  BERM.  BFRM*  ,F10  3) 

4700  WRITE(6,739)SFACE 

4800  739  FORMAT ( 2X .  ' THE  SLOPE  OF  THE  BEACH  FACE.  SF  ACE  *  . F 10  4) 

4900  WR  I  TE  (  6  ,  "T38  )  DIAM 

5000  738  FORMAT! 2X.  THE  SEDIMENT  DIAMETER.  DIAM*  .Flo  3) 

5100  WRITEI6.732) 

5200  780  F0RMAT(2X.  SUPPLY  MMAX  (  THE  NO  OF  GROINS)  AND  THEIR  I  -l.OC',/) 

5300  UCRIT* 16  3  *  SORT ( D I AM*0  00328 > 

5400  C*THE  NO  OF  MULTIPLE  GROINS. MMAX  MUST  BE  GIVEN  THEIR  X  LOCATIONS 
5500  READ! 5 . 779 )  MMAX 

5600  779  FORMAT! 13) 

5700  DO  760  MM,  MMAX 

5800  C*I JET  REPS  LESSER  I -VALUE  ADJACENT  TO  STRUCTURE 
5900  760  RE AD ( 5 , 779 )  IJET(M) 

6000  WRITE (6. 759)  ( M . I JE T ( M ) . M* 1 , MMAX ) 

6100  759  FORMAT ( 2X , ' THE  NUMBER  .15.'  GROIN  IS  IOCATED  AT  GRID  .15) 

6200  WR ITE(6,732) 

6300  C*CONVERT  TO  RADIANS 

6400  C *F I RST  MUST  GIVE  Y  COORS  POSITIONS  AND  DEPTHS 

6500  C'FIRST,  MUST  SET  UP  AIL  OF  THF  DEEP  VALUFS 
6600  WR !TE(6,773) 

6700  773  FORMAT! 2X, "NOW  ENTER  THE  VALUE  OF  ADEAN") 

6800  READ! 5 , 774  )ADE  AN 

6900  774  FORMAT! F 10  4) 

7000  WR I TE ( 6 , 749 )  ADc  AN 

7100  749  FORMAT ( 2X ,  1  THE  VjLUE  OF  ADEAN*  '  .  F  10  4.  IN  THE  EO  H*AV2/3  ) 

7200  WRI TE ( 6 . 732 ) 
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7  300 
7400 
7500 
7600 
7  700 
7800 
7900 
8000 
8100 
8200 
8300 
8400 
8500 
8600 
8700 
8800 
8900 
9000 
9100 
9200 
9300 
9400 
9500 
9600 
9700 
9800 
9900 
10000 
10100 
102 00 
10300 
10400 
10500 
10600 
10700 
10800 
10900 
1  1000 
1  1  ioo 
1  1200 
11300 
1  1400 
1  1500 
1  '600 
1  1700 
1  1800 
1  1900 
12000 
12  100 
12200 
12300 
12400 
12500 
12600 
12700 
12800 
12900 
13000 
13  100 
13200 
13300 
13400 
13500 
13600 
13700 
13800 
13900 
14000 
14100 
14200 
14300 
14400 


WR I TE ( 6 , 772 ) 

772  FORMAT ( 2X . "READ  IN  THE  SPACE  STE P . T I  ME  ST EP " . / ) 

READ! 5 , 775  I  OX.OELT 
775  FORMAT! 2<F 10  3 ) ) 

WR 1 TE ( 6 . 737 )  OX 

737  FORMAT! 2X .' THE  VALUE  OF  THE  LONGSHORE  SPACE-STEP.  DX  =  '  ,  F 10  3) 
WR I T  E ( 6 . 736 )  DELT 

736  FORMAT! 2X THE  TIME-STEP  IN  SECONDS.  DELT  =  ' . F 10  3) 

DATA  CHANGE/ 1..2..3..5  .7.11.14.17  .25  .32. 808 .10*0  0/ 

DO  220  J  = 1 . JMAX+2 
DO  220  1=1 . IMAX 
220  DEEP! I , J)=CHANGE( J) 

DATA(HC( I ) . 1  =  1 .8  )/1  87.0.5.0  35 .  25 .  2 1 .  20.  19.  19/ 

DATA! TC! I  )  . I  =  1 , 8)/2  . 3  . 4 ,  . 6  , 8  .  10  .12  .14  / 

DO  200  J= 1 . JMAX+2 
DO  200  I = 1 , IMAX 

200  Y( I ,  d+1 )  =  (0.5* (DEEP! I . J* 1 )+DEEP( I . J M/ADE AN ) *  *  1  5+  V ( I . 1 ) 

WR I TE ( 6 , 732  I 

r******t******»************* 

C*WE  WILL  ALWAYS  REOUIRE  Y ( I , JMAX+2 )  TO  COMPUTE  DY  AND  YBAR . 

C*WE  WILL  ALWAYS  REOUIRE  DEEP! I . JMAX+2 )  TO  COMP  SEDIMENT  TRANSPORT 

£«•**»*•*•*«***«********•** 


WRITE<6.734) 

734  FORMAT ( 2X .■ THE  BOUNDARY  Y-VALUES.  1=1. IMAX  ARE  AS  FOLLOWS’./) 

WR I TE ( 6 . 80 1 )  (Y(1.J),J=1. JMAX+2 ) 

WRITE (6 .801)  ( Y< IMAX . J) . J  =  t , JMAX  +  2 ) 

WR I TE ( 6 , 732 ) 

WRITE ( 6 . 735 ) 

735  FORMAT! / . 2X .  'THE  DEPTHS  8E  TWEEN  CONTOURS  ARE  AS  FOLLOWS  ./) 

WRITE (6. 801 )  (DEEP! 1 , J ) , d= 1 . JMAX+2 ) 

WR I TE ( 6 . 732  ) 

801  F0RMAT(2X, 10(F8  2)) 

DO  2  1=1. IMAX 

2  YZERO! I )*Y( I , 1 )-(BERM/SFACE ) 

C  +  W I L L  COMPUTE  THE  EOUIL  WIDTH  BETWEEN  CONTOURS.  HERE 
DO  1  1=1. IMAX 
WEO( I ,  1  )  =  Y( I  .  1 ) - Y ZERO ( I ) 

DO  1  d=  1  .  JMAX 
IMO  NE  II  GO  TO  32 
Y  TEMP  1 =0 . O 
GO  TO  33 

32  YTEMP1  =  ( (0  5*(DEEP( I  . J- 1 )+DEEP< I  . J)  )  J/ADEAN) •• 1  5 

33  YTEMP2=( (O  5* (DEEP! I , JI+OEEP! I , J+ 1 )) )/AOEAN)+* 1  5 
WEO! I , J+ 1 )=YTEMP2-Y TEMPI 

1  CONTINUE 

C*LET ' S  STORE  THE  ORIG  VALUES  TO  COMPUTE  VOL  CHANGES  OVER  CONTOURS . LATER 
DO  796  1  =  1 . I  MAX  + 1 
YZEROO! I  )=YZERO( I ) 

DO  796  J=1. JMAX+2 
796  YORIG! I . J)=Y( I , J) 

£********#•************♦**************•#*♦****♦****+♦**++♦***•++*=»***++ 
C+REAO  THE  DISK  FILE  AND  TRANSFORM  PARAMETERS  INTO  MY  UNITS 

C***********«***«***jfjfttfliltf'fitffffi  »»f  iff  Ml  («♦  +  +  **  +  +  •  +  ♦•  +  •  +  *•  +  *  + 

C+ALL  ADJUSTMENTS  TO  WAVE  ANGLE . HE IGHT . CELER I TY . GROUP  VEL.  WILL  BE  MADE 
C**HERE.  AND  THRUOUT  THE  REST  OF  THE  PROG.  THEY  WILL  BE  AS  IF  OCCURRED 
C+++AT  WDEPTH1 

798  READ! 5, 799. END= 1000)  HS.T.ALPWIS 

799  FORMAT! 10X . 3F6  1 ) 

NT  I MES= 1 

NCHECK  =NUNI V+NT IME  S 
HGEN=0  707 107+HS 
S I GMA  =  TWOP I /T 
G=  32  17 
CO=G*T/TWOPI 
ELO=CO*T 

IFIT.LE  2.0)  GO  TO  797 
HCC=0 . 23 
00  444  1=2.7 
T2  =  TC( I  ) 

IF!  T  GT  T2>  GO  TO  444 
T 1 =TC( I  -  1  ) 

DELTAT  =  T2-T  1 


67 


(4  500 
14600 
14700 
14800 
14900 
15000 
15100 
15200 
15300 
15400 
15500 
15600 
15700 
15800 
15900 
16000 
16  100 
16200 
16300 
16400 
16500 
16600 
16700 
16800 
16900 
17000 
17  100 
17  200 
17300 

1  7400 
17  500 
17600 

17  700 
17800 
17900 
18000 

18  100 
18200 
18300 
18400 
18500 
18600 
18700 
18800 
18900 
19000 
19  100 
19200 
19300 
19400 
19500 
19600 
19700 
19800 
19900 
20000 
20100 
20200 
20300 
20400 
20500 
20600 
20700 
20800 
20900 

2  1000 
2  1100 
2  1200 
2  1300 
2  1400 
2  1500 
2  1600 


OT  =  ( T-T  1  )/D£lTAT 
DTT  =  (  T2-T)/DELT 
HCC=HC ( I)*DT+HC(I-1)‘DTT 
GO  TO  446 
444  CONTINUE 
446  CONTINUE 

I  F  ( HGEN . LT . HCC )  GO  TO  797 

ANGGEN  =  A  L  PW I S  *  T WOP  1/360 

c. . . . . . . . 

CALL  WVNUM( WOE  PTH , T , DUMKK ) 

C  *  ANGGEN , HGEN , CGEN , CGGEN  REPRESENT  THE  WAVE  ANGL E . HE  I GHT  , C E L E R 1 T V  AND 
C  * 'GROUP  VEL ( RESPECT  )  OF  THE  SPECIFIED  WAVE  INPUT  AT  A  DEPTH.  WDEPTH 
CALL  WVNUM( 11.C.T, DUMKKK 1 
C 1 1 = TWOP I / ( T -OUMKKK ) 

CGI  1=0  5*C 1 1 *( 1 . +< 2 . ‘OUMKKK* 1 1  0/SINHI2  *DUMKKK*11  0))) 

CGEN  =  T WOP  I  /  ( T • DUMKK I 

CGGEN=0  5  *  CGEN  * (  1  ‘(2  ‘DUMKK ‘WDEPTH/S I NH( 2  ‘OUMKK'WDEPTH ) > ) 

CALL  TRANS 

797  IF1NCHECK  NE  NUNIV)  NUNIV=NCHECK 
709  GO  TO  798 
1000  CONTINUE 
STOP 
END 

c.  . . . . . . 

SUBROUTINE  TRANS 

C*THIS  SUBROUTINE  WILL  COMPUTE  SEO'IMENT  TRANSPORT 

COMMON/ A/  C<  60. 20 , RK ( 60 . 201 . Y ( 60. 20 ) .DEEP! 60. 20) . AL  PHAS1 60. 20  I 
COMMON/ AA/V ZERO (60) 

COMMON/BB/WEO( 60,20) 

COMMON/B/  THETA(60, 20) .QXTOT ( 60 ) ,  OLOANGI 60 . 20 ) ,  DVI60.20) 
COMMON/C/  H(60. 20) , CG< 60. 20) .HOLD ( 60, 20) . HB ( 60 . 20 ) . YB (  60 ) 

COMMON/N  USED/UUSE . T . CO , CGEN , CGGE N , ANGGEN . DX . BERM . THE T AO ( 10) . MM Ax 
COMMON/D/S  I  GMA.G.ELO,  JMAX,  I  MAX.  PI  .TWOP!  .  P  102.  HGEN.  I,JET<  101  .  SJ£TT> 
COMMON/ E /RHO. RHOS . POROS . CONST . TKSI 
COMMON/F/ADEAN, REPOSE , D I  AM 
COMMON /G/ IBREAK(60) .HNONBRt  20) 

C0MMON/P/HB0(6O) ,0EEPB(60) 

COMMON/ZZZ/NTIME 
COMMON/AA A/DEL T. NT IMES 
COMMON/COUNT /NUN  IV 

DIMENSION  YOLO (60. 20) . R ( 60 . 20 ) . S ( 60 , 20 ) . HC ( 60 . 20 ) . QY ( 60 . 20 ) . YD  I SS ( 
*  60.20) 

DIMENSION  RHS 1 ( 60 . 20 ) , S3 ( 60 . 20 ) . THE  TAB ( 60 . 20 ) . ANGLDC(60. 20 ) 
DIMENSION  0  I S  TR ( 60 , 20 ) , AWARE (60. 20 
C‘*‘*.‘*..‘‘.*.**‘..*.*‘****. 

£******.*.*.».*.*.. .*•**.*•***♦*.*..*.»*»*..***»**.*•»*.••• ............. 

c. ....... ......  NOTE  ■  SIZE  OF  ABAND  AND  XL  HAVE  TO  BE  CHANGED 

C*. ............  ACCORDING  TO  JMAX*1‘JMAX  AND  JMAX *  1 . RE SPE CT 

c. .............  CHANGE  REQ'D  AT  7040  AND  18650 

c. .........................  . . . . 

*  ) ,BMATRX( 432 ) , ABAND ( 432 . 19 ) . QX ( 60 , 20 ) . XL ( 432 , 10) . C0NST6 ( 60.20) 

COMMON/MP/  RKB ( 60 ) , HB I ( 60 ) . DE E PB I ( 60 ) 

COMMON/ EXPL/OY EXP (60. 20) , Y IMP (60. 20) 

DIMENSION  SANGLE ( 20 ) 

C'LET'S  ZERO-OUT  ALL  OF  THE  DIMENSIONED  MATRICES 
DO  1000  J= 1 , JMAX+2 
SANGLE ( J ) =0 . 0 
DO  1000  1=1, IMAX+2 
YOLD ( I . J)=0.0 
R( I , J ) =0  0 
S  (  I  , J ) =0 . 0 
HC ( I . J ) =0  -  0 
QY( I , J)=0  O 
YD  I  SS  (  I  ,  J  )  =0  0 
RHS  1  (  I , J ) =0  O 
S3 (  I  . 0 ) =0  0 
THE  T  AB ( I  . J)=0  0 
ANGLOC ( I , J ) =0  0 
0 1 S T R (  I  . J ) =0  0 
AWARE!  I  .  J  )  =0  O 
QX  (  I , J ) =0  0 
CONST 6 ( I  . J ) =0  0 
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2  1700 
2  1800 
2  1900 
22000 
22  100 
22200 
2  2  300 
2  2400 
2  2  500 
2  2600 
2  2  700 
22800 
22900 
2  3000 
2  3  100 
2  3200 
23300 
2  3400 
23500 
2  3600 
23700 
23800 
23900 
24000 
24  100 
24200 

24  300 
24400 
24500 
2  4  600 
2  4  700 
2  4  800 
2  4  900 
25000 

25  100 

25  200 
25300 
25400 
2  5500 
25600 
25700 
25800 
25900 
26000 

26  100 
26  200 
26300 
26400 
26500 
26600 
26700 
26800 
26900 
27000 
27  100 
2  7200 
27300 
27400 
2  7500 
27600 

27  700 
2  7800 
2  7900 
28000 

28  100 
28200 
28  300 
28400 
28500 
28600 
28700 
28800 
28900 


1000  CONTINUE 
RHO- 1  99 
RHOS -  5  14 
POROS  =  0  40 
CONS  T  =0  77 
CAPPA-O  78 
TAU-0  25 

TKS l  -  ( CONST'RHO'SQRT ( G  > )/( (RHOS -RHO ) ♦  (  1  0- POROS ) *  16  O' SORT  (  CAPPA ) 
C*  OX (1,0)  IS  the  transport  beiween  THE  (I.J+1)  AND  (I.J)  CONTOURS. 

C'THE  DO  1  LOOP  SIMULATES  TIME - T I  ME -DE L T 'NT  I  ME S 

COFF -0.00001 
GAMMA -RHO ' G 
DO  i  NT  I  ME  -  1  ,  NT  I  ME  S 
VJNIV'NUNIV' 1 

C ' 1  HE  MATRICES  ABANO  AND  BMATRx  MUST  BE  “ZEROED  OUT" 

K  =0 

00  26  I -2 .  IMAx  -  1 
DO  26  J- 1 . UMAX 
K-K*  1 

8MA  TRX(K)-0  0 
DO  26  L -  1 , JMA X  + 1 'UMAX 
26  ABANO ( K . L 1-0  O 

XNT I  ME  =  1 . 0'<NT IME ) 

CALL  PREOIF 

C  *  SMOOTHING  OF  THE  WAVE  ANGLE. THETA,  IS  RED  TO  ACCT  FOR  DIFF  EFFECTS. 
CALL  SMOOTH! THETA ,  IMAX , UMAX ,  I  JET  . S JETT  Y , MMAX . Y ) 

CALL  OTRAN 

C'FIRST  THE  LONGSHORE  SEDIMENT  TRANSPORT  WILL  BE  DISTRIBUTED 
C* ' ' 'ACROSS  THE  SURF  ZONE... 

CC= 1  25 

C"'QX(I,J)  WILL  BE  DETERMINED  By  SUBTRACTING  FROM  THE  INTEGRAL 
C  *  *  0  F  ox  r  ROM  OEEP(t,J-tl  TO  INFINITY,  THE  INTEGRAL  OF  OX  FROM  DEEP(I.J) 
C"*TO  INFINITY,  in  this  WAY  THE  SEDIMENT  TRANS  FROM  JMAX  out  gets 
C"  *  INCLUDED  IN  OX ( I  , JMAX I  TO  INCLUDE  THE  SWASH  TRANS.  WHEN  J- 1 
C'WE  WILL  SUBTRACI  FROM  2  TO  INFINITY  FROM  1.0 
C-LOOP  FOR  VALUES  WHICH  ARE  HELD  CONST  AND  STORED 
THETAB! 1 , 1 )=0  5’ ( THETA! 1 . 1 )»0  0) 

R( 1 . 1 )=0  5/(DX*(DEEP( 1 . 1 ) 'BERM/2  )) 

00  290  I -2 . IMAX 

R(I,1)=0  5/(DX*(DEEP( I . 1 ) 'BERM/ 2  )) 

C'  THETAB! I . 1 )=0  25* ( THETA! I  .  1 I'THETA!  I  -  1 .  1 )'0  ‘O  ) 

THE  TAB! I .  1 )=0.5'<  THETA! I ,  1 )* THETA! I- 1,  1 )) 

C'NO  NEED  TO  COMPUTE  PROP  ANGLE  AT  STRUCTS  BECAUSE  OX  -0  0  AT  IJETIMt'l 
ANGLOC! I,1)-ATAN((Y(I,1)-Y(I-1,1) )/0X ) 

C'HBO! I JET(M)*1 )  IS  PROPERLY  SET  IN  THE  SUBROUTINE  OTRAN 

0tSTR(I.1)*1.0-E*P(-((0EEP(I.1>"1.  5'HB0(  Il'ADEAN"  1  51/ 

'  (CC'DEEPB!  I  )"  1  5))"3) 

DISTR!  I  ,  1  ) -D I STR (  I  ,  1  )  *  TKS  I  'HBO(  I)  "2  5 
DO  290  J-2.JMAX 

R ! I . J ) =0  5/(DX' (DEEP! I . J) -DEEP! I , J- 1 ) ) ) 

THE  TAB (  I  . J ) -0  5*(THETA(  I  , Jl'THETA!  I  I.J)) 

ANGLOC! I , J ) -  A  T  AN ( (Y( I ,J)-Y(I-1,J))/0X) 

DISTR!  I  ,J)=EXP(  -  (  (DEEP!  I  ,J  -  1  )"1  5'HB0(  I  )' ADEAN"  1  5  )  /  (  CC  *DE  E  PE>  (  1  ) 
*  "1  5  )  )  "3  ) -EXP!  -  (  (DEEP!  I  ,  J)"  1 . 5'HBQ<  I  )*ADEAN"  1  ,  5  )  /  ( CC  ' 

'  DEEPB!  I  )"  1  5)  )  *  *  3  ) 

DISTR!  I , J ) -0  I S  TR!  I  .Jl'TKSI 'HBO! I )* '2  5 
290  CONTINUE 

00  301  J- 1 . JMAX 
DO  301  I -2 , IMAX 

aware(i.j)-delt'r(i,j)*(qx<i,j)-qx(i*i,j)+qy(i,j)-qy(i,j*i))*y(!  j 

'  ) 

si-2  'SIN! THE  TAB!  I  . J ) ) 'COS! THETAB!  I , J  )  ) * ( -  1  '2  '(COS! 

'  ANGLOC!  I  ,  J  )  )  )  "2  ) 

S2-C0S! 2  'THETAB! I  . J) ) 'COS ( ANGLOC ( I.J) ) / ( SORT ( DX  "  2* 

'  (  Y(  I  .  J)  -  Y(  I  -  1  ,  J  )  )  "2  )  ) 

S3! I , J)-S2'0rSTR( r . J) 

IF! SJETTY  EQ  o  O)  GO  TO  307 

DO  325  M= 1 .MMAX 

IF! I  NE  IJET(M)* 1 )  GO  TO  325 

IF(THETAO(M)  GE  00)  ISIDE-I JET(M) 

IF( THE  T  AO ( M )  LT  0  O)  I S I OE - I JE T ( M ) ' 1 

YSEA-0  5*( Y( I  SIDE . J) 'Y( I  SIDE , J' 1 ) ) 

Y SHORE -O  5'! Y(  ISIOE ,J)  +  Y( ISIOE ,J- 1 ) ) 
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2 9000 
29  '00 
29200 

2  9  300 
29400 
29900 
29900 
29700 
29800 
29900 
30000 
30100 
30200 
30300 
30400 
30500 
30600 
30700 
30800 
30900 

3  1000 
3  1100 
3  1200 
3  1  300 
3 '400 
3  1500 
3  1600 
3  1700 
3  1800 
3  1 900 
3  2000 
32  '00 
32200 
3  2  300 
32400 
32500 
3  2600 
32700 
32800 
32900 
3  3000 
33100 
33200 
3  3300 
33400 
33500 
33600 
33700 
3  3800 
33900 
34000 
34  100 
34200 
34  300 
34400 
34500 
34600 
34  700 
34800 
34900 
35000 
35100 
35200 
35300 
35400 
35500 
35600 
35700 
35800 
35900 
36000 
36  100 
36200 


I F (  Y  SE  A  G1  SJETTv  AND  YSHORE  GT  S  JE  TTY)  GO  TO  302 

I E ( ¥  S  E  A  GT  S JE  T  T  ¥  ANO  YSHORE  LE.SJETTY)  GO  TO  298 

C 'BECAUSE  A  NO  FlOW  B  C  IS  USED  ALONG  THE  STRUCT,  NO  ATTN  WAS  PAID 

C  * • T  0  GETTING  PROP  F  R  VALUES  CF  ANGLOC .  T HE T AB , D I STR . E TC 
S31 I , J  > -O  O 
OISTRI I ,Jl-0  O 
GO  TO  302 
325  CONTINUE 
GO  TO  302 

C  *  •  • A80VE .  ALL  PARAMETERS1I  E  . S 1 . S2 , S3 , THETAB . DISTR , ANGLOC ) 

C  *  * • ARE  COMPUTED  AS  IF  THE  STRUCT  IS  NOT  THERE  THE  B  C  AT  THE 
C  *  *  *  S  T  RUC  T  TIP  ASSUMES  QX  COMPUTEO  AS  IF  NO  STRUCT  PRESENT  AND  THEN 
C***BYPASSES  ACCORDING  TO  "RATIO" 

298  RATIO=(YSEA-SJETTY)/( Y SEA -YSHORE ) 

S  3 (  I  . JI-S3!  I  .JI-RATIO 
OISTRI I ,  J  )  =  D  I  S  T  R  I I .  J  )  •  P  A  T  I  0 
302  RHS 1(I.J)*DISTR<I.J)*S1-S3II.J)‘(Y(t.J)-Y(I~1,J)) 

301  CONTINUE 

CALL  BREAK! IMA* .JMAX ) 

C  *  TO  DETERMINE  DECAY  OF  C0NST6I I . J ) . NEED  WAVE  NO  AT  BREAKING. 

DO  754  1=1. IMAX* 1 

754  CALL  WVNUMI DEEPBI I  I  )  , T , RKB(  I  ) ) 

CAUSING  SHIELD’S  OIAG.Y  AXIS-0  05  8  ( T AUO-RHO-C *U *  *  2 ) . GET  UCRIT( FT/SEC ) 
UCRIT-16  3*SQRT(0IAM*  00328' 

00  750  I  - ' ,  IMAX* 1 
C0NST6I  I ,  1  1-COFF‘DX 
DO  750  J- 2. JMAX ‘2 

L • CONS  T6 (  I  ,  J )  GOES  W/  QYII.J)  WHICH  IS  ASSOC  W/  DEEP(I.J-I) 
IF(OEEPII.J-I)  LE  DEEPBIII))  GO  TO  751 
C-HFRE  must  CAUSE  COFF  TO  DECAY  (WE’RE  BEYOND  SURF  ZONE) 

UMAXB-HBI I  I  ) ‘G‘T*RKB( I)/( 2 . ‘TWOPI *C0SH(RK8( I ) -DEEPBI ( I ) ) ) 

UMAX -H(  !.J-1)‘G*T‘RK<I.J-1)/(2  *  TWOP I *COSH( RK(  I.J-1)*DEEP(I.J-1))) 
IFIUCRIT  LT  UMAX  ANO  UCRIT.LT  UMAXBI  GO  TO  749 
CONS T 6 (  I , J ) =0  0 
GO  TO  750 

749  TOP-0  01*H<  I  .  J-  1  )“3‘SIGMA“3/(SINH!RK(  I  ,0-  1  )*OEEP(I  .d-  1  ))“3) 

SOT -DEEP!  I  . J ■  1  I  • I  0  625* TWOPI *G“ 1  5*0  78 • • 2  * ADE AN* *  1  5  + 

*  I  0  01 ‘HBI I  I  ) • *3 'SIGMA* *3/1 DEEPBI ( I )*( SINHIRKBI I ) ‘DEEPBI ( I ) ) )  “  3 )  I ) 
CONS  T  6 1  l , j)=Dx*C0FF‘T0P/B0T 

GO  TO  750 

75'  CONSTCI I ,J1=C0FF*DX 

750  CONTINUE 
K  =  0 

C**PUT  INTO  BANDED  FORM  USING  THE  ALGORITHM  A ( M . N ) - >B ! M , NN )  WHERE 
C“*NN  =  KB*1  M  +  NI KB  IS  THE  NUMBER  OF  LOWER  COO  I AGONALS! - JMAX . HERE ) ) 

DO  304  1=2. IMAX -  1 
00  304  J  = 1 , JMAX 
K  =  K*  1 

AWARE(I.J)-AWAREII.J) *DELT ‘RHS 1 ( I , U ) *R( I . J ) -DELT*R< I . J ) ‘RHS 1 1  I  *  1 . J 

*  )‘DEL  T*R( I  . J ) ‘C0NST6I I . J)‘WEO( I.U)-DELT*R(I . J)*C0NST6( I  . J*  1  ) * 

*  WFOI  I . J ‘  1  ) 

YDUM-YZEROI I ) 

IF! J  NE  1 )  YDUM=Y( I , J- 1 ) 

AWARE ( I . J (-AWARE ( I . J  )  *0E L T ‘ R (  I , J ) ‘C0NST6 ( I , J ) *0 . 5‘ ( YDUM- Y ( I . J)) 

*  -OELT'RI I . J) ‘C0NST6I I . d* 1 ) *0  5‘ ( Y ( I , J ) - Y ( I , J+ 1 ) ) 

U=DELT*R( I . J)‘S3( I . J) 

V-OEL  T‘R(  I , J  )  »  S3 1  I  +  1 , J ) 

71-0ELT*R(I . J)*C0NST6( I . J ) *0 . 5 
Z2-DELT‘R( I . J ) *C0NST6( I ,d+ 1 )‘C  5 
C ‘NOW  WILL  SET  UP  THE  MATRICES  ABANO  AND  BMATRX 
ABANO! K , JMAX* 1  ) -  1  . 0+U*V*Z 1 *Z2 
IF(  I  NE  2  )  GO  TO  305 
AWARE! I . J) -AWARE  I  I . J)+U‘Y( 1-1 ,d) 

GO  TO  310 

305  ABANDIK.  1  ) -  - U 

310  I F ( I  NE  IMAX-1)  GO  TO  306 

AWARE! I  . J) -AWARE! I  , J ) * V* Y ( IMAX . J ) 

GO  TO  311 

306  ABANDIK , JMAX* 1* JMAX )= -V 
3  1  1  I  F  (  J  NE  1  )  GO  TO  307 

ABANDIK , JMAX* 1  )  -ABANDI K , JMAX* 1 ) - Z 1 
AWARE!  I  .  1  1-AWARF I  I  .  1  )*Z1*|YZER0I  I  ) - Y 1  I  .  1  )  ) 

GO  TO  312 
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36300 

307 

ABANO ( K . UMAX ) - -Z 1 

l 

1 

36400 

312 

IF  (  J  NE  JMAX  )  GO  TO  308 

36500 

AWARE! I . J)=AWARE( I . J )*Z2*Y< I . JMAX* 1 ) 

1  • 

1 

36600 

GO  TO  309 

\ 

36700 

308 

ABANO ( K . JMAX  +  2 )  =  -  Z2 

1 

i 

36800 

309 

BMATRX(K)=AWARE( I , J) 

! 

36900 

304 

CONTINUE 

3  7000 

KMAX  =K 

37  100 

C  * ‘CALL  I MSL  ROUTINE  LEQT1B  TO  SOLVE  THE  BANDED  MATRIX. 

37200 

CALL  LEQT IB ( ABAND . KMAX , JMAX . JMAX , 432 . BMATRX , 1 . 432 . 0 . XL . I ER ) 

37  300 

C-NOW 

GIVE  YS  THEIR  NEW  VALUES  STORING  OLD  VALUES  IN  YOLD 

37400 

K  =  0 

; 

37500 

DO  315  1=2. IMAX- 1 

i 

37600 

YOLD! I  . JMAX* 1  >  =  Y(  I  . JMAX* 1 ) 

! 

37700 

DO  315  J=  1 . JMAX 

37800 

K  =  k  *  1 

37900 

YOLD! I . J!=Y( I . J) 

' 

38000 

Y ( I , J ) =BMATRX ( K 1 

38  100 

315 

CONT INUE 

38200 

DO  320  J= 1 , JMAX* 3 

* 

38300 

YOLD! 1 , J)«Yl 1. J) 

i/ 

38400 

320 

YOLD! IMAX . J)=Y( IMAX , J) 

38500 

C*WILL  USE  ABBOTT'S  DISSIPATIVE  INTERFACE  TO  RID  HIGH  FREO  OSCILLATIONS 

38600 

DO  650  J= 1 , JMAX 

38700 

00  650  I =2 . IMAX- 1 

*  •' 

38800 

YD  I SS ( I  . J)  =  TAU*Y( I-1,J)  +  (1.-2.*TAU)*Y(I,J)*TAU*Y(I*1.J) 

38900 

IF! SJETTY . EQ.O.O)  GO  TO  650 

39000 

DO  649  M= 1 , MMAX 

! 

39100 

I F ( I  NE . I JET(M)  AND . I . NE  I JET(M)+ 1 )  GO  TO  649 

39200 

IF! Y( I JET(M) . J)  GT  SJETTY .OR  Y( I JET(M)* 1 . J) .  GT . SJETTY )G0  TO  649 

39300 

IF( I . EQ. IJET(M) JYOISS! I . J )  =  TAU* Y ( I  -  1 . J )  +  ( 1 . - TAU ) * Y (I . J ) 

39400 

IF!  I  EQ  ( I JET(M)* 1  )  IYOISS! I  . J)=TAU*Y( I* 1  . J)*(  1  . -TAU)*Y( I  . J ) 

39500 

649 

CONTINUE 

39600 

650 

CONTINUE 

39700 

DO  651  J= 1 , JMAX 

39800 

DO  651  1=2. IMAX- 1 

39900 

651 

Y! I . J)=YDISS( I , J) 

40000 

C-THIS  LOOP  WILL  STORE  THE  IMPLICIT  V  VALUES  REQ'D  TO  COMP  QY8QX 

40100 

DO  40  1=1 . IMAX+1 

• 

40200 

DO  40  J= 1 , JMA X * 3 

40300 

40 

Y IMP!  I  .  J  ) =Y( I . J  ) 

40400 

C*THIS  LOOP  WILL  EXPLICITLY  MOVE  CONTOURS  SEAWARD  IF  REPOSE  EXCEEDED 

40500 

KOUNT  =  0 

40600 

SL0PEM=TAN(0  9-REPOSE) 

40700 

DO  48  1=1 , IMAX 

40800 

43 

KOUNT =KOUNT* 1 

40900 

IF (KOUNT  GT  50000)  GO  TO  4 1 

4  lOOO 

C  *  LET 

US  COMPUTE  ALL  THE  SLOPE S ( PSLOP )  FOR  EACH  CHANGE  IN  DEPTH 

4  1100 

DO  47  J= 1 , JMAX* 1 

4  1200 

DUM= -BERM/2  0 

4  1300 

I F ( J  NE  1)  DUM=DEEP( I , J- 1 ) 

4  1400 

DELH  =  0. 5* (DEEP! I  . J* 1 ) *DE EP ( I , J ) ) -0  5  * ( DE EP ( I . J)  +  DUM) 

4  1500 

PSLOP=DELH/( Y( I . J* 1 ) - Y ( I . J ) ) 

4  1600 

47 

S ANGLE ( J ) = AT  AN! PSLOP) 

4  1700 

C * F I NO  THE  MIN  NEG  SLOPE  ANGLE  OR  THEN  THE  POS  SLOPE>REPOSE  OR  FORGET  IT 

' 

4  1800 

ASLOPM= -  1  0E50 

4  1900 

ASLOPP  =0 . 0 

42000 

DO  46  J= 1 . JMAX* 1 

42  100 

I F ( SANGLE ( J )  GT  00)  GO  TO  45 

42200 

I F ( SANGLE ( J )  GT  ASLOPM) A SLOPM= SANGLE ( J ) 

4  2300 

IF( ASLOPM. EQ  SANGLE! J))  JM=J 

42400 

GO  TO  46 

42500 

45 

IF(SANGLE(J)GT  REPOSE  AND  SANGLE(J)  GT  ASLOPP )ASLOPP=SANGLE ( J ) 

42600 

I F ( ASLOPP  EQ  SANGLE! J) )  UP=J 

1 

42700 

46 

CONTINUE 

42800 

IF(ASLOPM. EO. -  1  0E50  AND  ASLOPP  EQ.O.O)  GO  TO  42 

42900 

I F ( ASLOPM  EQ.-1.0E50)  GO  TO  44 

| 

4  3000 

DUM* -BERM/ 2 . 

43100 

IF(JMNE.I)  DUM'DEEP! I . JM- 1 ) 

4  3200 

ALTER=( (0  5/SL0PEM-(0EEP( I , JM* 1 ) -DUM) ) - ( Y ( I , JM* 1 ) - Y ( I . JM ) ) ) / 

43300 

( 1  0*( (DEEP! I . JM*1 )-DEEP( I , JM) )/(DEEP( I , JM)-DUM) ) ) 

4  3400 

Y ( I . JM* 1 )=Y( I ,JM* 1 J+ALTER 

71 

i  | 

1 

43500 
4  3600 
43700 
4  3800 
4  3900 
44000 
44  100 
44200 
44300 
44400 
44500 
44600 
44700 
44800 
44900 
45000 
45100 
45200 
45300 
45400 
45500 
45600 
45700 
45800 
45900 
46000 
46100 
46200 
46300 
46400 
46500 
46600 
46700 
46800 
46900 
47000 
4  7  100 
47200 
4  7300 
47400 
47500 
47600 
47700 
4  7800 
47900 
48000 
48100 
48200 
48300 
48400 
48500 
48600 
48700 
48800 
48900 
49000 
49100 
49200 
49300 
49400 
49500 
49600 
49700 
49800 
49900 
50000 
50100 
50200 
50300 
50400 
50500 
50600 
50700 


Y ( I . JM ) * Y ( I . JM)-( ALTER- (DEEP ( 1 . JM  + 1 ) - DE EP ( I . JM ) )/ ( DEEP( I . JMJ-DUM) ) 
QYEXP( I .  JM  + 1 )=OYEXP( I , JM  + 1 ) +DX/DE LT ‘ALTER* < DEEP ( I . JM* 1 ) -DEEP( I , JM ) 
*  ) 

GO  TO  43 
44  CONTINUE 

DUM* -BERM/2 . 

IF(JP.NE.I)  DUM=DEEP( 1 , JP- 1 ) 

ALTER=((0  5/ SLOP EM* ( D£EP( I , JP+ 1 )-DUM) )-( Y( I , JP+1 )-Y( I , JP ) ) )/ 

*  ( 1  0+( (OEEP< I . JP* 1 )-DEEP( I , JP) )/(DEEP( I . JP ) -OUM ) ) ) 

Y ( I , JP* 1 )=Y( I . JP* 1 )+ALTER 

Y( I . JP)=Y( I . JP)-( ALTER *( DEEP ( I , JP* 1 )-DEEP< I . JP ) )/ ( DEEP ( I . JP)-DUM) ) 
OYEXP( I , JP+1 )=QY EXP ( I , JP+1 )+DX/DELT ‘ALTER* (DEEP( I . JP* 1 ) -DE EP ( I . JP ) 

*  ) 

GO  TO  43 

42  WEQ( I . JMAX+ 1 ) =Y ( I . JMAX+ 1 ) - Y( I , JMAX ) 

48  CONTINUE 

C*IF  WE  GET  SENT  HERE.  LOOP  444  WILL  CATCH  THE  CROSSED  CONTOURS 
41  CONTINUE 

C*NOW  WE  CAN  COMPUTE  OX'S  AND  OY ' S 1 
DO  318  1=2.1  MAX 

C*ALL  IMPLIC  AND  EXPLIC  MOVEMENT  OF  Y ZERO  WILL  BE  TAKEN  CARE  OF  HERE 
0Y( I . 1 )  =  -BERM*OX  * ( Y ( I .  1 )-YOLD( 1,1) )/DELT 
YZERO( I )=YZERO( I )+( Y(I . 1 )-YOLD( 1,1)) 

319  DO  318  J=1 , JMAX 

QX( I , J)=RHS1 ( I ,J)-S3( I,J)*YI MP ( I ,J)  +  S3( I ,J)*Y IMP( I  -  1 , J ) 

318  OY ( I . J+1 )=C0NST6( I . J+ 1 ) * ( 0 . 5* ( Y IMP ( I . J)+YOLD( I . J ) - Y IMP ( I . J+ 1 ) 

*  -YOLD( I , J+1 ) )+WEQ( I , J+1 ) ) 

DO  323  J=1 .JMAX 
QX( 1 . J ) =QX ( 2 , J ) 

323  0X( IMAX+1 . J)=OX( IMAX, J) 

C'TOTAL  QYS  WILL  BE  COMP  FROM  IMPLIC  AND  EXPLIC  VALUES. THEN  ZERO  QYEXP 
DO  39  1=1 . IMAX+1 
DO  39  J= 1 . JMAX+3 
QY< I . J)*QY(I . J)+QYEXP< I , J) 

39  QY EXP (  I  , J ) =0  O 

C*THIS  CHECK  WILL  BOMB  THINGS  OUT  IF  CONTOURS  HAVE  CROSSED. 

DO  444  11=1, IMAX 
DO  444  JJ=1 . JMAX 

C *  I F  CONTOURS  CROSS  AT  ANY  TIME  WANT  PROGRAM  TO  STOP  1 
I F ( Y ( I  I , JJ) .LT . Y( II , JU+1 ) )  GO  TO  444 
WR I TE ( 6 , 103) 

WRITEI6.*/)  NUN IV 
DO  150  J=1 , JMAX 

150  WRITEI6, 100)  (0X( I , J) , I = 1 . IMAX ) 

DO  151  J= 1 , JMAX 

151  WR ITE(6, 101 )  (OY(I.J). 1=1, IMAX) 

DO  152  J=1 , JMAX 

152  WRITE(6, 1 OO )  ( Y( I , J ) , I = 1 . IMAX ) 

103  FORMAT (2X, 'THE  CONTOURS  HAVE  CROSSED  AND  SOMETHING  IS  WRONG'./) 

DO  19  J=1 .JMAX 

19  WRITE ( 6 , lOO )  ( YOLD( I . J) , I = 1 . IMAX) 

GO  TO  445 
444  CONTINUE 

WR ITE(6,*/)  NUNIV 

C*THE  FOLLOWING  STATEMENT  DETERMINES  AT  WHAT  FREO  EVERYTHING  IS  WRITTEN' 
I F ( MOD ( NUN IV, lO)  .  NE .0)  GO  TO  1 

C*LET ' S  WRITE  ALL  OF  IT  OUT. 

WR I TE ( 6 , 926 )  NUNIV 

926  F0RMAT(2X, 'THE  TOTAL  ELAPSED  NUMBER  OF  TIME-STEPS.  NUNIV*  ’,15,/) 

800  FORMAT ( 2X, 1 4 ( F8  4  )  ) 

C*  DO  900  1=1, IMAX 

C*900  WR I TE ( 6 . 800 )  ( THETA( I . J ) , J= 1 . JMAX ) 

C*  DO  903  J-1.JMAX+1 

C  *903  WR I TE ( 6 , 80 1 )  DEEP(I.J) 

C*  DO  906  1=1, IMAX 

C*906  WR I TE ( 6 . 800 )  ( H( I , J ) . J* 1 , JMAX ) 

C*  DO  755  J* 1 .JMAX 

C*755  WRI T£(6 , 800 )  ( C0NST6 ( I , J ) . I = 1 , IMAX ) 

801  FORMAT (2X, 14(F82)) 

WRI TE ( 6 . 107) 

107  F0RMAT(/,2X, 'THE  LONGSHORE  TRANSPORTS . OX ,  FOLLOW) 

DO  15  J= 1 , JMAX 

15  WRITE! 6 , 10O )  (OX(I ,J) , 1=1 , IMAX ) 
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BA>Oo  MBOI  I  I  *C  AA*A*A  *0A  A**ft  (  I  1 

BAAOO  I  Ofttgi  I  t  ANO  PA  |A*ftl  t  A  ft*  1 1  A  BA  (OMfolAO  AiiltBOINo  AO  AmA  BAVA  OIB 

BBAOO  (••  AA  tMf  At  Bill'  tllBA  III*  DlfAAO 

BAftlftl  |  A  |  ft.lA  t  A  a  Ag  O  "A  (MA  l)A  I 

ftAAOO  00  A  M* I  MMAa 

BAAOO  I A  |  |  NA  |  *l>  A  |  M 1  •  I  A  IMA  10  A 

ftftBOO  I  ••AMA  ABANAAOBA  |AA(t  BAVA  ft  Bill  BA  (OMPlUAO  UMAi  A"A  BAVA  All  I'BOI-  fttl'l 
BAAAim  |  f  |  AM|  A  AOI  M )  OA  O  01  (Ml  Dl  IA 

ft  Aland  (lAffBI  I  l*(M(  |,IA  A  A  M  I  •  A  ,1*  1  I*  (1  All*  A  l.'A  A|MD  I  .11  ll"n  a  iAAI'AI"i'  ft 

BAAOO  I  BBA  Aft  A  |  I  -  ABBA  AMI  I  ,'A  A  ( B  I  •  |  ) 

ftltiai  'MA  10  IJ 

BAAOO  II  OA  I  A*B  A  I  A  *  A  Alt  I  .'A  •  A  M  1  ,M  I  A  •  IM  A  A*  A  l.'A  A  A  M I  ...  '  I  «  «0  }ft  I  AI*A>A  I  •  D*  ft 

BAAOO  IBBAAKA  I  I  -  | BBA  Aft  I  |  ,iA  A  (Ml  I 

BABOO  I  >  'Btgil  t -OAAfBI  I  ID  AM*A 

BAftOO  (MA  AO  I 

ftAAOO  ft  (ONAINOA 

BABOO  00  All  I 

BAftOO  >  (ONAINOA 

BBOOO  A  I  ONI  IAAAiA 


58100 
58200 
58300 
58400 
58500 
58600 
58700 
58800 
58900 
59000 
59100 
59200 
59300 
59400 
59500 
59600 
59700 
59800 
59900 
60000 
60100 
60200 
60300 
60400 
60500 
60600 
60700 
60800 
60900 
61000 
61  100 
6  1200 
6  1300 
6  1400 
61500 
61600 
61700 
61800 
6  1900 
62000 
62100 
62200 
6  2300 
62400 
62500 
62600 
62700 
62800 
62900 
63000 

63  100 
6  3200 
63300 
63400 
63500 
63600 
63700 
63800 
63900 
64000 

64  100 
64200 
64  300 
64400 
64500 
64600 
64700 
64800 
64900 
65000 
65100 
65200 


C  *  I F  THE  OFFSHORE  WAVE  HT  IS  ZERO.  NEVER  GET  TO  HERE 

C'HOWEVER  IF  THE  H  IS  SUCH  THAT  IT  WOULD  BREAK  INSHORE  OF  Y ( I . 2 ) 

C*0EEP8( I  )  WOULD  STILL  BE  ZERO  AND  OISTR(I.J)  WOULD  BLOW-UP 
DO  20  1=1. IMAX 

IF(DEEPB( I )  GT.OO)  GO  TO  20 

OE  EPB ( I )  =  (H( 1 .  1)‘DEEP(I,1)“0. 25/CAPPA ) * *0 . 8 

HBQ( I )=CAPPA*DEEPB! I ) 

20  CONTINUE 

HBQ( 1 )=H8Q(2 ) 

H8Q( IMAX  + 1 ) =HBO( IMAX ) 

DEEPB(  1 ) =DE  EPB ( 2 ) 

OEEPBI IMAX+1 )  =OEEPB ( IMAX ) 

RETURN 

END 

. . . . . . . . . . . 

SUBROUTINE  BREAK! IMAX , JMAX ) 

C‘ROUTINE  WILL  DETERMINE  HB  AND  DEEPB  ON  THE  GRID  LINES  RATHER 
C*  THAN  BETWEEN  THEM  REQ'O  FOR  COFF  BEYOND  SURF  ZONE 

COMMON/ A/  C(60. 20) ,RK(60, 20) . Y ( 60 , 20) , DE EP( 60 . 20 ) . ALPHAS ( 60 . 20 ) 
COMMON/C/  H(60. 20) .CG(60, 20) ,HOLD( 60 . 20 ) , HB( 60 . 20 ) . YB ( 60 ) 
COMMON/MP/  RKB ( 60 ) , HB I ( 60 ) . DE EPB I ( 60 ) 

CAPPA=0 . 78 
DO  1  1=2. IMAX 
DO  2  JJ=1 .UMAX 
U= JMAX- JU+ 1 

IF(H( I , J) , LT  HB! I ,J) )  GO  TO  2 

DEEPBI! I )=( (H( I . J+ 1 )*DEEP( I . J+ 1 ) * *0  25 )/CAPPA ) • ‘O  8 
HBI ( I )=CAPP A ‘DEEPBI ( I  ) 

C  “  ‘ONCE  THE  HEIGHT  &  DEPTH  AT  BREAKING  ARE  FOUND.  GO  TO  NEXT  GRID-LINE 
GO  TO  1 
2  CONTINUE 
1  CONTINUE 

DO  20  1=1. IMAX 

IF(DEEPBI( I )  GT .0.0)  GO  TO  20 

DEEPBI!  I  )  =  (H(  I  .  1  (‘DEEP!  I  .  1  >  “0 . 25/CAPPA  )  “0  8 

HB I ( I )*CAPPA‘DEEPBI ( I ) 

20  CONTINUE 

DEEPBI! 1 )=0EEPBI(2) 

DEEPBI! IMAX-*- 1 )=0£EPBI( IMAX ) 

HBI( 1 )=HBI(2) 

HBI ( IMAX*1 ) =HB I ( IMAX ) 

RETURN 

ENO 

C*. ..**.*.*•.*.*.*.**...***•****.*.*.***.**♦*...*..»......*.***. ........ 

SUBROUTINE  REFRAC ( UBEGIN , UENO . NPTS , IBEGIN.  I  END , 1START.M) 

COMMON/ A/  C!60.20) , RK( 60 . 20) . Y ( 60 , 20 ) . DEEP ( 60 . 20 ) . ALPHAS ( 60 . 20 > 
COMMON/ A A/ Y ZERO (60) 

COMMON/B/  THETA ( 60 . 20 ) . QXTOT ( 60 ) .  OLDANG! 60 . 20 ) ,  DY(60,20) 
COMMON/C/  H(60. 20) ,CG(60, 20) .HOLD (60, 20) . HB ( 60. 20 ) . YB ( 60 ) 

COMMON/N  USED/ JUSE  ,  T  ,  CO  ,  CGEN ,  CGGEN  ,  At  IGGEN  ,  DX  .  BERM  ,  THET  AO!  10)  .MMAX 
COMMON/D/S IGMA.G.ELO, UMAX, I  MAX, PI . TW )P I .PI02.HGEN. IUET(  10) .SUETTY 
COMMON/ G/IBREAK! 60 ) ,HNONBR(20) 

COMMON/ZZZ/NT I  ME 

DIMENSION  UBEGIN! 60 ) . UEND! 60 ) 

c... .............  THIS  SUBROUTINE  WILL  DETERMINE  H  AND 

c. ...............  THETA  AT  THE  MID  PT  OF  Y  VALUES 

C  “  • T  AU  IS  THE  FACTOR  WHICH  RECOUPLES  THE  REFRACTION  EOS. SEE  ABBOTT 
T AU  =  0  25 

C ‘MUST  PRESCRIBE  THE  WAVE  ANGLE  AT  THE  OUTERMOSTCONTOUR  BOX 
C‘SNELL'S  LAW  WILL  BE  USED  TO  START  THINGS  OFF 

C ‘THET A ! I , U )  WILL  BE  AT  AREA'S  CENTER  AND  WILL  USE  Y(I.U)  IN  NEG  Y-DIR 
C‘WILL  INITIALIZE  ALL  THETA'S  USING  SNELL'S  LAW. 

DO  206  I ‘IBEGIN. I  END 

C‘ INITIALIZE  TWO  U-VALUES  BEYOND  UMAX. IF  IN  REGION  1. 

I F( UENO! I ) . EO ■ UMAX )  UINIT=2 

IF(UEND( I )  NE  UMAX )  UINIT-0 
DO  206  U= UBEGIN! I ) , UENO! I )+UINIT 
C ‘MUST  CORRECT  FOR  THE  CONTOUR  ORIENTATION.  ALPHAS. 

IF( I  NE . IBEGIN)  GO  TO  960 

ALPHAS!  I  .  U  )  ‘AT  AN(  (0.5*(Y(I*1,U)‘Y(I*1.U*1))0  5  *  (  Y  (  I..  U  ) 

»  ♦Yd.u+on/ox) 

GO  TO  962 
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65300  960  IF (I  NE  I  END)  GO  TO  961 

65400  ALPHAS! I , d ) = AT AN< ( 0  5* ( Y ( I  . 0 )  +  Y ( I . J+ 1 ) ) -0  5* ( V ( I  -  1 , J ) 

65500  *  ♦Y(I-1.d+1)))/DX) 

65600  GO  TO  962 

65700  961  ALPHAS!  I  .  d)*ATAN(  (O  5*<  Y(  1  +  1  ,d)+Y(  1+  1  .  d  +  1  )  )-0  5* 

65800  *  ( V(  I- 1 ,d)*Y< I- 1 ,d+ 1 ) ) )/<2 . *DX ) ) 

65900  962  DALPHA* ANGGEN- ALPHAS ( I . J ) 

66000  THE  T  A ( I , J  ) = ARS IN! ( C ( I , d  ) /CGEN ) *S IN! DALPHA ) > 

66100  C'MUST  GET  THETA  WRT  THE  X-AXIS. 

66200  THETA! I . J )  =  THE  T  A ( I ,d)  +  ALPHAS< I . J) 

66300  206  CONTINUE 

66400  C*NOW,  WE  MUST  COMP  THE  80UN  WAVE  HTS  SO  THE  HTS  CAN  BE  COMPUTED 

66500  C*WILL  USE  THE  EO.  •***♦*  DEL  DOT  (E»CG)=0.0 

66600  C*NOW  WE  WILL  CORRECT  THE  HT  FOR  SHOALING  AND  REFRACTION  TO  THE  B  C 
66700  C-WILL  ALSO  INITIALIZE  H'S  WITH  THESE  EQUATIONS  FOR  ENTIRE  ARRAY 
66800  DO  500  I  =  IBEGIN , I  END 

66900  CMNITIALIZE  TWO  d-VALUES  BEYOND  UMAX  IF  IN  REGION  1. 

67000  I F ( UEND ( I )  EO . UMAX )  JINIT=2 

67100  IF! JEND! I  )  .  NE  JMAX )  dINIT*0 

67200  DO  500  d  = JBEGIN! I ) . JENO( I ) +UINI T 

67300  H ( I , d)=HGEN*SORT(CGGEN/CG! I . d ) ) *  SORT < COS < ANGGEN ) /COS! THE T A ! I . 

67400  *  d ) ) ) 

67500  I F ( HB ( I . d ) . L  T . H ( I . d ) )  H( I , d ) =HB< I  , d ) 

67600  500  CONTINUE 

67700  C* - - - - - - - 

67800  (;**■•*******♦*♦**’••*****♦*♦**♦**»**•♦•*♦*»♦***♦*•****♦*♦****•*♦•♦*♦»♦••••♦ 

67900  C*L£T ' S  FILL  THE  DY  ARRAY. 

68000  C*DY  WILL  BE  INDEXED  AS  THE  THETA  TO  WHICH  WE  ARE  GOING. 

68100  DO  209  I  * IBEGIN . I  END 

68200  DO  209  d=dBEGIN( I ) *1 , dEND! I ) 

68300  DYlI ,d-1)=0.5'(Y(I.d-l)*Y(!,d) |-0.5*(Y( I ,d)*YI  I ,d+1 ) I 

68400  209  CONTINUE 

68500  NITERS* lOO 

68600  DO  lOO  NITER=1 .NITERS 

68700  SUMANG=0  O 

68800  C*DO  "60  LOOP"  GOES  FROM  2  TO  IMAX  IF  ISTART  * IBEGIN 

68900  C "DO  "60  LOOP"  GOES  FROM  IMAX-1  TO  1  IF  I S  T  AR  T= I END 

69000  DO  60  1 1  =  IBEGIN , I  END 

69100  C'MUST  HAVE  IT  SET  UP  SO  THAT  THE  KNOWN  BOUNDARIES  ANGLES  AREN'T  RECOMP 
69200  IF! ISTART  EO. IBEGIN)  I*!I 

69300  IF(ISTART.EQ  IBEGIN  .AND.  I. EO. IBEGIN)  GO  TO  60 

69400  IF! ISTART. EO. IEND)  I  =  I  END- 1 1  +  IBEGIN 

69500  IF( ISTART  EO. IEND  AND.  I.EQ.IEND)  GO  TO  60 

69600  C*ADX  EQUALS  ACTUAL  OELTA  X  ACROSS  SPACE  STEP 

69700  C "ONLY  ON  BOUNDARIES  WHERE  FORWARD  OR  BACKWARD  DIFFERENCING 

69800  IF! I  NE . IBEGIN)  GO  TO  6 

69900  ADX=DX 

70000  IP=I+1 

70100  IM*  I 

70200  GO  TO  12 

70300  6  IF! I  NE . IEND)  GO  TO  10 

70400  ADX  =  DX 

70500  IP*  I 

70600  I M* I  -  1 

70700  GO  TO  12 

70800  10  ADX=2  0*DX 

70900  IP*I+1 

7  1000  I M= I  -  1 

71100  12  CONTINUE 

71200  DO  40  d*dBEGIN( I ) ,dENO( I )- 1 

71300  C*WILL  GO  FROM  (dMAX-1)  TO  1  BECAUSE  THAT'S  THE  DIR  WAVE  COMES  IN  FROM 
71400  dd*dEND(I)-1-d+dBEGIN!I) 

71500  OLDANG! I ,dd)*THETA( I ,dd) 

71600  C*LOCATE  MIDPOINT  BETWEEN  TWO  ADdACENT  BLOCK  CENTERS 

71700  C ‘BECAUSE  THETA'S  dd-VALUE  IS  THE  SAME  AS  THE  FIRST  SHOREWARD  Y  VALUE 

71800  C*MUST  USE  dd,  dd+ 1 .  AND  dd+2  TO  COMPUTE  YBAR . 

7 1 900  YBAR*0  25* ( Y( I , dd ) *2 . 0* Y ! I , dd+ 1 )+ Y ( I . dd*2 ) ) 

72000  C*LOCATE  APPROPRIATE  INDICES  ON  IP  AND  IM  GRID  LINES 
72100  IMINUS* -  1 

72200  I  PLUS*  *■  1 

72300  CALL  LOC ( IM, dd , dO IM . dSIM . YBAR , IM] NUS ) 

72400  CALL  LOC ( I P , dd . dOI P . US  I P . YBAR , I  PLUS ) 

72500  C*NOW  USE  THE  CONSERVATION  OF  WAVES  EQUATION . 
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72600 
72700 
72800 
72900 
73000 
73100 
73200 
7  3300 
73400 
73500 
73600 
73700 
73800 
73900 
74000 
74100 
74200 
74300 
74400 
74500 
74600 
74700 
74800 
74900 
75000 
75100 
75200 
75300 
75400 
75500 
75600 
75700 
75800 
75900 
76000 
76100 
76200 
76300 
76400 
76500 
76600 
76700 
76800 
76900 
77000 
77100 
77200 
77300 
77400 
77500 
77600 
77700 
77800 
77900 
78000 
78100 
78200 
78300 
78400 
78500 
78600 
78700 
78800 
78900 
79000 
79100 
79200 
79300 
79400 
79500 
79600 
79700 
79800 


PART 1C*RK( I ,  dd+1 )*SIN(THETA< I , dd+ 1 ) ) 

PART2*-0Y(I .dd)/A0X 

C*WILL  LINEARLY  INTERPOLATE  TO  DETERMINE  RK*COS( THETA )  AT  1+1  AND  1-1. 
C*IF  NO  ADd  SHOREWARD  PT  EXISTS,  PUT  IN  ZERO  FOR  TERMS  IN  GOV.  EO. 

IF ( dSIM . NE .0)  GO  TO  301 

PART38*0 . 0 
GO  TO  302 

301  TOPIM»RK( IM.dOIM- 1 ) *COS( THET  A{ IM, dOIM- 1 ) ) 

BOTIM*RK( IM.dSIM)*COS(THETA(IM.dSIM) ) 

TOT ALB =0  5*(Y(IM,UOIM)+Y(IM,dOIM- 1 ) ) -0  5* ( Y( IM, JSIM+ 1 )+Y( IM, JSIM) ) 

DUMB*0.5*(Y(IM.dOIM)+Y( IM.dOIM- 1 )  )  - YBAR 

PART3B* ( ( TOT ALB -DUMB ) * ( TOP IM-BOT IM  J/TOTALB ) +BOT IM 

302  IF(dSIP.NE .0)  GO  TO  303 
PART3A*0 . 0 

GO  TO  304 

303  TOPIP-RK( IP.dOIP- 1 )»COS(THETA< IP.dOIP- 1 ) ) 

BOTIP*RK( IP,dSIP)*CQS( THET  A( IP.dSIP)) 

TOTALA-O. 5* (Y( IP.dOIP )+Y( IP.dOIP- 1 ) ) -O  5«< Y ( IP , dSIP* 1 )+Y ( IP , dS I P ) ) 
DUMA*0.5*(Y( IP.dOIP )+Y( IP.dOIP- 1 ) ) - YBAR 
PART3A*((T0TALA-0UMA)*( TOPIP-BOTIP  ) /TOTAL A)+BOTIP 

304  PART3-PART3A-PART3B 

C*NOW  MUST  FIND  RK*SIN( THETA )  FOR  1+1  AND  1-1  AT  d+ 1 
YBARP-0.25»(Y( I . dd+ 1 )+2 . *Y( I ,dd+2 )+Y( I ,dd+3) ) 

CALL  LOCUM,  dd+1  .dPOIM.dPSIM.YBARP,  IMINUS) 

CALL  LOC ( IP, dd+1 , dPOIP . dPSIP. YBARP , I  PLUS) 

I F ( dPSIM . NE . 0 )  GO  TO  305 
PART 1B*0 . O 
GO  TO  306 

305  TOPM*RK( IM.dPOIM- 1 ) *SIN( THETA( IM,dP0IM-1 ) ) 

BOTM-RK( IM, dPSIM) *SIN( THETA ( IM.dPSIM) ) 

T0TB*O. 5+( Y( IM.dPOIM )+Y( IM.dPOIM- 1 ) ) -0 . 5* ( Y ( IM , dPSI M+ 1 )+ 

•  YUM, dPSIM)) 

DUMPB*0.5*( Y (IM.dPOIM )+Y( IM.dPOIM- 1 ) ) -YBARP 
PART  IB* ( ( TOTB -OUMPB ) * ( TOPM-BOTM )/TOTB )  +  BOTM 

306  IF ( dPSIP  NE.O)  GO  TO  307 
PART  1 A  *0 . O 

GO  TO  308 

307  TOPP*RK( IP.dPOIP-1 )*SIN( THETA ( IP,dP0IP-1 ) ) 

B0TP*RK(I P.dPSlP)+SIN( THET A  UP, dPSIP ) ) 

TOTA.O.5*(Y(IP,dP0IP)+Y(IP,dPOIP-1 ))-0.5*(Y( I P . dPS IP+ 1 ) +Y ( IP. dPSIP 

*  )) 

DUMP A  *0 . 5  * ( Y ( IP.dPOIP)  +  Y( IP.dPOIP- 1 ) ) -YBARP 
PART  1  A* ( ( TOTA-DUMPA ) * ( TOPP-BOTP )/TOTA )+BOTP 

308  PART 1-TAU*PART1B+( 1.-2. • T AU ) *PART 1C+TAU+PART 1 A 
IF(dPSIM.EQ.0)PART1*( 1 . -TAU) *PART 1C+TAU+PART 1A 
I F ( dPSI P . EO.OJPART 1*TAU+PART 1B+( 1 . -TAU)*PART1C 
ARG* ( ( PART 1+PART2*PART3)/RK( I . dd ) ) 

C*IF  THE  ROUTINE  IS  TO  BLOWUP, USE  SNELLS  LAW 
I F ( ABS( ARG )  LE  1  0 )  GO  TO  41 
ARG* ( C ( I , dd )/C( I , dd+ 1))*SIN( THETA ( I . dd+ 1 ) ) 

IF(ARG.GT  1.0)  ARG* 1 . 0 
THETA ( I ,dd)*ARSIN( ARG) 

GO  TO  42 

41  THET  A(I,dd)*ARSIN(APG) 

42  THETA ( I , dd ) *0 . 5* ( THETA ( I , dd ) +OLDANG( I , dd ) ) 
SUMANG*SUMANG+(ABS(THETA( I ,dd)-OLDANG( I . dd ) ) ) 

40  CONTINUE 
60  CONTINUE 

C*MUST  EdECT  IF  WE  HAVE  REACHED  AN  ACCEPTABLE  ITERATION  ERROR 

C*IF  THE  SUM  OF  THE  ABSOLUTE  VALUE  OF  ANGLE  CHANGES  DURING  AN  ITERATION 

C*  AVERAGES  LESS  THAN  0.02  DEGREES  PER  GRID  ITS  CLOSE  ENOUGH. 

I F ( SUMANG . LT . ( NPTS+O . 0035 ) )  GO  TO  215 
IF(NITER  GE .50)  GO  TO  215 
100  CONTINUE 

WRI TE (6 , 803 ) 

215  CONTINUE 

C* ITERATION  LOOP  FOR  THE  WAVE  HEIGHT. 

00  501  NITER*1 .NITERS 
SUMH*0  0 

00  510  1 1  * IBEGIN, I  END 

C ‘MUST  HAVE  IT  SET  UP  SO  THAT  THE  KNOWN  BOUNDARIES  HTS .  AREN'T  RECOMP 
IFUSTART  EQ  IBEGIN)  I  - 1 1 

I F ( ISTART . EO. IBEGIN  AND.  I. EO  IBEGIN)  GO  TO  510 
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79900 
80000 
80100 
80200 
80300 
80400 
80500 
80600 
80700 
80800 
80900 
81000 
81  100 
81200 
81300 
8  1400 
81500 
81600 
81700 
81800 
81900 
82000 
82100 
82200 
82300 
82400 
82500 
82600 
82700 
82800 
82900 
83000 
83100 
83200 
83300 
83400 
83500 
83600 
83700 
83800 
83900 
84000 
84  100 
84200 
84300 
84400 
84500 
84600 
84700 
84800 
84900 
85000 
85100 
85200 
85300 
85400 
85500 
85600 
85700 
85800 
85900 
86000 
86100 
86200 
86300 
86400 
86500 
86600 
86700 
86800 
86900 
87000 


I F ( I  ST ART. EQ  I  END)  I ■ I  END- 1 1 +  IBEGIN 
IF( ISTART . EQ. IEND  AND  l.EQ.IEND)  GO  TO  510 
C * ADX  EQUALS  ACTUAL  DELTA  X  ACROSS  SPACE  STEP. 

C‘ONLY  ON  BOUNDARIES  WHERE  FORWARD  OR  BACKWARD  DIFFERENCING. 

IF( I  NE . IBEGIN)  GO  TO  503 
ADX-DX 
IP‘I+1 
I  M=  I 

GO  TO  505 

503  I F  (  I  .  NE  .  I  END  )  GO  TO  504 

AOXsOX 

IP=I 
IM=I-1 
GO  TO  505 

504  ADX  *2 . 0*DX 
IP‘I+1 

IM= I  -  1 

505  CONTINUE 

DO  502  d=dBEGIN( I ) , dENO( I )- 1 
dd=dENO(I ) -  1 -d+dBEGIN( I ) 

HOLD( I . JJ)-H( I . Jd) 

YBAR=0.25‘(Y( I ,dd)+2 ,0‘Y( I , dd+ 1 )  +  Y ( I , dd+2 ) ) 

CALL  LOC< IM.  dd.  dOIM.  dSIM,  YBAR  , IMINUS) 

CALL  LOC(IP.dd,dOIP.dSIP.YBA«, I  PLUS) 

PART 13=(H( I . dd+1  )“2 . ) *CGf I .dd+1 )*COS( THETA ( I , dd* 1 ) ) 

PART2=DY( I , dd )/ AOX 
IF(dSIM.NE.O)  GO  TO  311 
PART  4B=0 . 0 
GO  TO  312 

31 1  TOPIMH=(H( IM.dOIM-  1  )  “2  . )*CG( IM.dOIM- 1 )*( SIN( THETA ( IM , dOIM- 1 ) ) ) 

BOT IMH= ( H( IM.dSIM)**2 . )*CG( IM . dSIM ) ‘SIN( THETA ( IM , dSIM ) ) 

T0TALB=0  5*( Y( IM, dO IM ) *Y ( IM, dOIM- 1 ) ) -0 . 5‘ ( Y ( IM . dSIM* 1 )  +  Y ( IM . dSIM ) ) 
DUMB=0.5‘(Y( IM. dOIM )*Y( IM.dOIM- 1 ) ) -YBAR 

PART4B*((T0T  ALB -DUMB ) * ( TOPIMH-BOT 1MH )/TOTALB )*BOTIMH 

312  IF(dSIP.NE.O)  GO  TO  313 
PART4A -0 . O 

GO  TO  314 

313  TOPIPH=(H(IP,dOIP- 1 )**2 . ) »CG( I P . dOI P- 1 ) *SIN( THET A( I P , dOI P- 1 ) ) 
BOTIPH=(H( IP.dSIP>*«2  ) *CG( IP . dSI P ) *S IN( THET A( IP.dSIP) ) 
TQTALA=0.5‘(Y( IP . dOIP )*Y < IP , dOIP- 1 ) ) -0 . 5» ( Y ( IP , dSIP* 1 )*Y( IP . dSIP ) ) 
0UMA=0  5*(Y(IP. JOIP)*Y( IP.dOJP-1 ))-YBAR 

PART4A  =  ( (TOTALA -DUMA )  • ( TOP  I  PH- BOT IPH)/TOTALA )+BOT I PH 

314  PART4=PART4A-PART4B 

YBARP=0 . 25* ( Y ( I , dd* 1 ) *2  *  Y ( I ,dd+2 )  +  Y ( I ,dd*3) ) 

CALL  LOC(IM.dd*1.dPOIM,dPSIM.YBARP, IMINUS) 

CALL  LOC( IP. dd*1 . dPOIP.dPSIP . YBARP , I  PLUS) 

IF(dPSIM  NE  O)  GO  TO  315 
PART  12=0.0 
GO  TO  316 

315  TOPMH=(H( IM.dPOIM- 1 )**2 )*CG( IM, dPOIM- 1 ) *COS( THETA ( IM . JPOIM- 1 ) ) 
BOTMH* (M(IM,dP5IM)**2)*CG{ IM , dPS IM ) *COS ( THETA ( IM. dPS IM ) ) 

TOTB*  5*(Y( IM,dPOIM)*Y( IM.dPOIM- 1 ) )- . 5 • ( Y ( IM . dPS IM* 1 )+ Y ( IM.dPSIM) ) 
DUMPB=0 . 5* ( Y ( IM. dPOIM)+Y( IM.dPOIM- 1 ) ) -YBARP 
PART 12»( ( TOTB -OUMPB ) * ( TOPMH-BOTMH )/TOTB )+BOTMH 

316  I F ( UPS  IP  NE  0 )  GO  TO  317 
PART  11=0  0 

GO  TO  318 

317  TOPPH* ( H( IP.dPOIP- 1 )**2)*CG( IP.dPOIP- 1 )*COS(THETA( IP.dPOIP- 1 ) ) 
BOTPH* (H(IP,dPSIP)**2 ) *CG( IP . dPSIP)*COS<  THETA ( IP . dPS IP ) ) 

TOT A-  5*(Y( IP.dPOIP )*Y( IP.dPOIP- 1 ) )-  5 * ( Y ( I P . dPS I P* 1 )*Y( IP . dPSIP ) ) 
DUMP A =0  5*< Y( IP ,dPO IP )+Y( IP.dPOIP- 1 ) ) -YBARP 
PARTI W (TOTA-OUMPA)*( TOPPH-BOTPH)/TOTA )+BOTPH 

318  PART 1H*TAU*PART 12*( 1.-2. *T AU ) ‘PART  1 3*T AU*PART 1 1 
IF(dPSIM  EQ  OfPART 1H*<  I  -TAU)*PART 13*TAU‘PART 1 1 
IF ( dPSIP. EQ. OPART 1H‘TAU*PART12*<  1 . - T AU ) ‘PART  1 3 

ARG‘( (PART 1H*PART2*PART4 )/(CG( I . dd ) ‘COS( THE TA( I , dd ) ) ) ) 

C* I F  THERE  IS  TO  BE  AN  INVALID  SORT, USE  LINEAR  SHOALING. 

I F ( ARG . GE  0. )  GO  TO  44 

ARG* ( CG( I , dd* 1 ) ‘COS( THET  A( I . dd* 1 ) ) )/( CG< I . dd ) ‘COS( THET A ( I .dd) ) ) 
IF(ARG.LT  0.0)  ARG‘0.0 

H( I ,dd)«H( I ,dd* 1 ) ‘SORT ( ARG ) 

GO  TO  45 
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87  100 
87200 
8  7  300 
87400 
87500 
87600 
87700 

8  7800 
87900 
88000 

88  100 
88200 
88300 
88400 
88500 
88600 
88700 
88800 
88900 
89000 
89100 
89200 
89300 
89400 
89500 
89600 
89700 
89800 
89900 
90000 
90100 
90200 
90300 
90400 
90500 
90600 
90700 
90800 
90900 
91000 

91  100 

9  1200 
91300 
91400 
91500 
91600 
91700 
91800 
91900 
92000 

92  100 
92200 

92  300 
92400 
92500 
92600 
92700 
92800 
92900 
93000 

93  100 
93200 
93300 
93400 
93500 
93600 
93700 
93800 
93900 
94000 

94  100 
94200 
94  300 


44  H( I , dd)*SQRT(ARG) 

45  H( 1 , JJ)’0  5*(H( 1 . dd)+H0LD< I . dd) ) 

HNONBR (  dd ) =H( I , dd ) 

C* IBREAKd )  =  dd.  THEREFORE  dd  WILL  BE  LEEWARD  SIDE  OF  GRID  AT  INIT  BREAK 
IF(HB(I.dd)  .LT.  H(I,dd)  .AND.  HB (I . dd+ 1 ). GE . HNONBR < dd+ 1 ) ) 

*  IBREAK( I ) =dd 

IF(HB< I ,dd)  LT  H( I .dd) )  H< I ,dd)=HB( I . dd ) 

SUMH  =  SUMH  +  ABS ( H ( I ,dd)-HOLD( I . dd ) ) 

502  CONTINUE 
510  CONTINUE 

I  BREAK! IEND)=IBREAK( I  END- 1 ) 

I  BREAK! I  BEG  I N )  * IBRE AK( IBEGIN+ 1 ) 

IF(SUMH.LT . (NPTS*0. 01 ) )  GO  TO  507 
IF(NITER.GE  50)  GO  TO  507 
501  CONTINUE 

WRITEI6.803) 

507  CONTINUE 

802  F0RMAT(2X,4(F 15.5),////) 

803  F0RMAT(2X. "AFTER  NITERS  ITERATIONS,  CONVERGENCE  WAS  NOT  REACHED") 

804  FORMAT (2X, "THE  WAVE  HT .  ROUTINE  CONVERGED  IN.  NITER*  ".15.//) 

805  F0RMAT(2X. "THIS  IS  MY  CHECKING  WRITE  STATEMENT") 

806  F0RMAT(2X , "THF  WAVE  ANGLE  ROUTINE  CONVERGED  IN.  NITER*  ".15.//) 
RETURN 

END 


C*********************************************************************** 


SUBROUTINE  DIFF ! RHONO . THETAO, ANGLE . AMP ) 

C****0IFFRACTI0N  ABOUT  SEMI  INFINITE  BREAKWATER  I PENNE Y -PR I CE ) 
P 1  =  3  14159265 

A6SS*SIN(0  5* (ANGLE -THETAO) ) 

ABSP=SIN(0  5*( ANGLE +THET AO) ) 

ABC =COS( ANGLE -THETAO) 

ABC 1*C0S< ANGLE* THETAO) 

XX*RHOND*ABC 

XXC=C0S(XX) 

XXS=SIN( XX ) 

XX 1 *RHOND*ABC 1 
XXC1=C0S(XX1 ) 

XXS1=SIN(XX1 ) 

AL=SQRT(RHOND/PI ) 

SIG=2  0*AL*ABSS 

S I GP* -  2 . 0* AL  * ABSP 

CALL  FRES(SIG,C.S.FR,FI ) 

CALL  FRES(SIGP.CP.SP.FRP.FIP) 

SUM 1=XXC*FR+XXS*FI +  XXC 1 +FRP+XXS 1  *  F I P 
SUM2*XXC*FI -XXS+FR+XXC 1 *F I P -XXS 1 *FRP 
AMP=SQRT(SUM1 **2+SUM2**2) 

RETURN 

END 


SUBROUTINE  FRES ( A . C , S . FR , F I ) 

C'FRESNEL  INTEGRAL  SUBROUT INE **• *AF TER  ABROMOWITZ  AND  STEGUN. 
Z*ABS( A) 

P02  *  1  5707963 

FZ*( 1 .0*0. 926*Z)/(2 .0+1. 792»Z+3 . 104*Z*Z) 

GZ  = 1 .0/(2. 0+4.  142  +  Z  +  3 . 492*Z*Z  +  6  670*Z*Z*Z) 

XX  *P02  +  Z  *  Z 
CZ'COS(XX) 

SZ=SIN( XX ) 

C*0  5-GZ*CZ+FZ*SZ 

S=0  5-FZ*CZ-GZ*SZ 

IF! A  GT.O.O)  GO  TO  50 

C=-C 


S*-S 

50  FR*0  5*( 1  O+C+S) 
FI=-0  5*(S-C) 
RETURN 
END 


SUBROUTINE  PREDIF 

COMMON/ A/  C(60. 20) ,RK(60, 20) . Y(60. 20) .DEEP (60. 20) . ALPHAS ( 60 . 20 ) 
COMMON/ AA/Y ZERO! 60) 

COMMON/B/  THETA(60. 20) ,0XT0T(60) .  OLDANG! 60. 20) .  DYI60.20) 
COMMON/C/  H(60. 20) ,CG(60, 20) .HOLD! 60. 20) . HB ( 60 . 20 ) . YB ( 60 ) 
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94400 
94500 
94600 
94  700 
94800 
94900 
95000 
95100 
95200 
95300 
95400 
95500 
95600 
95700 
95800 
95900 
96000 
96100 
96200 
96300 
96400 
96500 
96600 
96700 
96800 
96900 
97000 
97  100 
97200 
97300 
97400 
97500 
97600 
97700 
97800 
97900 
98000 
98100 
98200 
98300 
98400 
98500 
98600 
98700 
98800 
98900 
99000 
99100 
99200 
99300 
99400 
99500 
99600 
99700 
99800 
99900 
tooooo 
100100 
100200 
100300 
100400 
100500 
100600 
100700 
100800 
100900 
101000 
101 100 
101200 
101300 
101400 
101500 
101600 


COMMON/N  USED/JUSE  .  T  .  CO  ,  CGEN  .  CGGEN  .  ANGGEN  ,  DX",  BERM  ,  THE  T  AO  (  10)  .  MMAX 
COMMON/D/ S I GMA ,G . ELO , JMAX . I  MAX . PI , TWOPl .  PI02 . HGEN. I  J£T(  10  I . SJET  T  v 
COMMON/G/ I BRE  AK  <  60  > , HNONBR { 20 ) 

DIMENSION  J1(60) ,d2(60) . J1R£F(60) . J3REF(60) 

C'THIS  SUB  CALCS  WHERE  0 1  EE  RACE  I  ON  GOVERNS  AND  WHERE  REFRACT  GOVERNS 
C*IT  WILL  CALL  REERAC  FOR  OFFSHORE  ARE  A( OF  F  TIP  OF  STRUCTURE) 

C  *  THEN  IT  WILL  00  THE  SHAOOW  ZONE  USING  0 1 F  F ( IF  THETAO  NE  O  O) 

C*  IT  WILL  THEN  FINISH  THE  OTHERS  USING  REFRAC  AGAIN. 

C'LET'S  ZERO-OUT  THE  DIMENSIONED  ARRAYS 
DO  lOOO  1=1 . 1MAX+2 
J 1 ( I ) =0  O 
J2 ( I ) =0 . O 
J 1  RE  F ( 11=00 
1000  d3REF (  I  )  =0  0 

C'NOW.  LETS  FIND  C.CG.RK.HB,  AND  WVNUM 
DO  202  I  =  1  .  I  MAX 
DO  202  J  = 1 , JMAX-*- 2 
DEPTH=DEEP( I . d ) 

CALL  WVNUM ( DEPTH . T , DUMK ) 

RK ( I . d ) =DUMK 

C( I , J)=CO'TANH(RK< I . d ) *DE EP (  I  . d ) ) 

EN  =  0  5* (  1  0+( (2  *RK( I  ,d)‘OEEP( I . J) )/SINH(2 , *RK( I . J)*DEEP<  I  .d)  )  ) ) 

CG ( I , J)=EN*C( I ,d) 

H8 ( I , J ) *  0 . 78*DEEP( 1.0) 

202  CONTINUE 

C  *WI LL  ATTR  FB  AN  EQUAL  REACH  TO  EACH  SIDE  OF  EACH  M-GROIN 
DO  200  M= 1 ,MMAX 
10UML=1 

IF(M  NE  1)  IDUML=( I  JET < M )  + 1  JET ( M- 1 ) )/2 
I DUMR  = I  MAX 

IF(M  NE  MMAX )  I DUMR  =  ( I dET ( M )  + 1  JET < M+ 1 ) ) /2 
NPTS=0 

DO  1  I = IDUML , IDUMR 
DO  2  J«1.JMAx 

I F ( Y ( I . J  )  LT  SdETTY)  GO  TO  14 
d1(I  )  =  d 
J2( I ) =  JMAX 
GO  TO  1 5 

14  CONTINUE 
2  CONTINUE 

15  CONTINUE 

CMF  NO  STRUCT  IS  PRESENT* SJETTy=0  O) .  DO  REFRAC  THRUOUT  GRID  SYSTEM 
IF(SdETTY  EO.O  01  J 1  (  I  )  =  1 
1  CONTINUE 

DO  16  I  *  I DUML  .  I  DUMR 

C*  'REFRAC'  STARTS  ON  THE  NEXT  TO  LAST  J-CONTOUR . NOT  THE  LAST) 

00  16  J  =  Jl  (  I  ) . d2< I  ) -  1 

16  NPTS=NPTS+1 

C’WILL  NOW  DO  THE  REFRACT  FOR  THE  REGION  1  AREA 
C'ISTART  REPRESENTS  THE  DIRECTION  THE  SWEEPS  WILL  BEGIN  FROM 
C  * W I LL  USE  DUMMY  IMAX . I  JET . IdET* I  IN  call  SITS  SO  IBEGIN.IEND.  and 
C** • ISTART  WON'T  CHANGE  THEM  MUST  RESET  AFTER  EACH  CALL  REFRAC 
IMAXTdDUMR 
IJETT=IJET(M> 

IJETP1=IJET(M>* 1 
IOUMLL  =  I DUML 

IF(ANGGEN.GE  00)  CALL  RE FRAC ( d 1 . J2 , NPT S . I DUML l . IMAX T . I DUMLL . M ) 

I F  (  ANGGEN  LT  0  0)  CALL  RE FR AC ( d 1 , J2 . NPTS . I DUML L . I MAXT , I  MAX T . M ) 

IMAXT  =  IDUMR 

IdETT=IdET(M) 

IdETP1=IdET(M)+1 
IDUMLLdDUML 
JDUMN=d1( IJET(M) ) 

J0UMS=d1 ( IdETl M)+ 1 ) 

XDISTN=( IJET(M)- 1  0 ) *0X  FDX/2 

ELTIP*T*0.5*(C( IdET(M) . JDUMN)*C( I JE  T ( M ) ♦ 1 , JDUMS ) ) 

C*NOW  MUST  CHECK  THE  ANCLE  AT  THE  STRUCTURE'S  TIP  TO  SEE  WHERE  SHAD  ZONE 
C  *  I F  NO  STRUCT  PRESENT(SJETTY=0.0) •  FJTHER  REFRAC/DIFF  UNNECESSARY 
I F ( SdETTY  EQ  0  0)  GO  TO  13 

THE  T  A0( M ) =0  5*( THETA( IdET (M) , JDUMN ) *  THE T A ( IdET(M)< 1 .JDUMS ) ) 

HINC=0  5*(H( IJET(M), JDUMN )*H( IdET (M)* 1 , JDUMS) ) 

I F ( THETAO(M) ) 10.  H . 12 

C*THIS  SECTION  HANOLES  RFFR4C/0IFF  IF  THCTA0<0  0 
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101700 
101800 
101900 
102000 
102100 
102200 
102  300 
102  400 
102500 
102600 
102700 
102800 
102900 
103000 
103 100 
103200 
103300 
103400 
103500 
103600 
103700 
103800 
103900 
104000 
104100 
104200 
104300 

104400 

104500 
104600 
1047QO 
104800 
104900 
105000 
105100 
105200 
105300 
105400 
105500 
105600 
105700 
105800 
105900 
106000 
106100 
106200 
106300 
106400 
106500 
106600 
106700 
106800 
106900 
10 '000 
107100 
107200 
107300 
107400 
107500 
107600 
107700 
107800 
107900 
108000 
108  100 
108200 
108300 
108400 
108500 
108600 
108700 
108800 
108900 


10  CONTINUE 

C-FIRST  ALL  OF  REGION  2  WILL  GET  REFRACTED 
NPTS=0 

DO  lOO  I  =  IJET(M»*1 .  I  DUMP 
J2( I  )  =  0 1 ( I  ) 

100  J 1 ( I  )  -  1 

DO  101  1=1 JET(M)* 1 . IOUMR 
DO  101  U  =  U1( I ) , J2 ( I  )  -  1 

101  NPT  S  =  NPT  S  + 1 
IMAXT* IOUMR 
IDUMLL* IDUML 
IJETT=IJET(M) 

IUETP1=IUET(M)*1 

CALL  REF  RAC ( J 1 . J2 ,NPTS.  IJETP1 , IMAXT .  IMAXT ,M) 

I  MAX  T  =  I DUMR 
I  UE  T  T  =  I  UE  T  (  M  ) 

IUETP1=IUET(M)* 1 
IDUMLL=IDUML 

C  *  NOW  MUST  00  REGION  3  OF  NEG  THE  T  AO  CASE-SHADOW  ZONE. 

DO  102  I  =  I DUML . I UE  T ( M  ) 

U2) I  )  =  U 1 <  I  ) 

102  J 1 ( I >  = 1 

00  103  I  =  I OUML . IUET(M) 

J  1  RE F ( I  )  =  1 

DO  104  U  =  U  1  ( I) . U2< I )+ 1 

XCOOR  =  ( I  -  1  . O) *0X 

YC00R=O  5* ( Y( I . U )  +  Y( I . U+ 1 >  ) 

ANGLE -AT  AN( ( XDI STN- XCOOR )/( SUET TY-YCOOR ) ) 

I F ( YCOOR  GT  SUETTY)  ANGLE =PI * ANGLE 
C-IF  MOST  SHOREWARD  PT  OUT  OF  SHAD  ZONE.  SO  ARE  THE  OTHERS  FOR  THAT  I 
IF(ABS(ANGLE )  GT  ABS( THETAO(M) ) )  GO  TO  105 
RAD * SORT ( (XDISTN-XC00R)**2+( SUETTY -YCOOR)* *2) 

RHOND*RAD*  TWOP I /ELT I P 
C  *0 1 FFRACT I  ON  TREATS  THE  POS  THETAO  CASE. 

THE*ABS(THETAO(M) ) 

CALL  01 FF ( RHONO . THE . ANGLE . AMP ) 

H( I . U ) =AMP*HINC 

angrao*-angle 

C  *WI LL  NOW  REFRACT  OIFF  WAVES  IN  THE  SHAD  ZONE  USING  SNELL'S. 
CTIP*ELTIP/T 

ALPHAS) I ,U)=ATAN( (0  5* ( Y < l +1 , J )+Y ( I + 1 . U+1 ) ) -0 . 5* 

*  ( Y( I  -  1 , U)*Y< I  -  1 ,u+ 1 ) ) )/( 2 . *DX ) ) 

I F ( I  EQ  IUET(M) (ALPHAS) I . U ) = ATAN) ( 0  5* ( Y < I . U )  +  Y ( I . U* 1 ) ) -O . 5 * ( Y ( I  -  1 

*  . U)*Y( I  -  1 ,U* 1 ) ) )/DX ) 

DAL  PHA - ANGR AD -ALPHAS) I  . U ) 

THETA) I . U)=ARSIN( (C( I . Ul/CTIP )*SIN(DALPHA ) ) 

THETA( I , U ) *  THE  T  A ( I ,U)*ALPHAS( I , U ) 

C • MUS T  CHECK  TO  SEE  IF  WAVE  WOULD  HAVE  BROKEN 

IF(HB)I.U)  LE  H< I . U )  AND  HB( I , U+ 1 )  GT  H( I ,U* 1 ) )IBREAK( I)*U 
IF(HBII.U)  LT  H(l.U))  H(I.U)-HB)I.U) 

104  CONTINUE 
GO  TO  103 

105  U 1REF ( I ) =U 

103  CONTINUE 

C  *  NOW  MUST  DO  REFRACTION  FOR  REGION  4 
NPTS*0 

DO  106  I = IDUML . IUET(M) 

DO  106  U  *  U  1  R  E  F (  I  ).U2(I )- 1 

106  NP  T  S  =  NP  T  S  +  1 
IDUMLl * IOUML 
IMAXT -  I DUMR 
IUETT-IUET(M) 

IUETP1*IUET(M)* 1 

CALL  REFRAC) U1REF ,U2 . NP T S . IDUMLL . IUETT . IDUMLL ,M) 

I DUMLI.  *  I  OUML 
I MAXT  *  I DUMR 
IUETT=IUET(M) 

IUETP1*IUET(M)* 1 
GO  TO  13 

C'THIS  HANDLES  REFRAC/DIFF  IF  THETAO  IS  00 
C  *  F  OR  THIS  CASE.  ONLY  THREE  REGIONS  EXIST 

11  CONTINUE 
NP  T  S*0 
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109000 
109100 
109200 
109300 
109400 
109500 
109600 
109700 
109800 
109900 
1 10000 
1 10100 
1 10200 
1 10300 
1 10400 
1 10500 
1 10600 
1  10700 
1 10800 
1  10900 
1  1  lOOO 
1 1 1 100 
1 1 1200 
1  1  1300 
1 1 1400 
1 1 1 500 
1 1 1600 
1 1 1 700 
1 1 1800 
1 1 1 900 
1  12000 
1 12100 
1 12200 
1 12300 
112400 
1  12500 
1 12600 
112700 
1 12800 
112900 
1  13000 
1 13100 
1 13200 
1 13300 
1 13400 
1 13500 
1 1 3600 
1 13700 
1 13800 
1 1 3900 
1 1 4000 
1 14100 
1 14200 
1 14300 
1 14400 
1 14500 
1 14600 
1 14700 
1 14800 
1 14900 
1 15000 
1 15100 
1 15200 
1 15300 
1 15400 
1 15500 
1 15600 
1 15700 
1 15800 
1 15900 
1 16000 
1  16  100 
1 16200 


DO  120  I = I DUMl , 1JEIIMI 
J2( I ) = J 1 ( I  ) 

120  J  1  (  I  )  =  1 

DO  121  I = IDUML . I OET( M) 

DO  121  J  =  J1(  I) . J2<  I  )  1 

121  NPTS=NPTS+1 
I  MAX  T  = 1 DUMR 

I DUML  L - 1 DUML 
I  JE  T  T  =  I  JE  T  (  M  > 

I OE TP  1  =  I JE T ( M  )  + 1 

CALL  REE RAC ( J1 . J2.NPTS. I DUML L. 1JETT . IDUMLL.M) 

IMAXT  =  IDUMR 
1J£TT=I JET(M) 

IJETP1=IJET(M)*1 
I DUML  L  =  I DUML 

DO  122  I = I  JE  T ( M )  ♦ 1 , IDUMR 
02 ( I  )  =  J1< I  ) 

122  J  1  (  I  )  =  1 
NPT  S  =  0 

DO  123  I  =  IJET  <  M ) *  1 ,  IDUMR 
DO  123  J  =  J1  < I  )  . J2< I) -  1 

123  NPTS=NPTS+ 1 
IMAXT  =  IDUMR 
IDUMLL  =  I DUML 
IJETT=IJET(M) 

IJETP1=IUET(M}+ 1 

CALL  REFRAC(J1.J2,NPTS.IJETP1.IMAXT,IMAXT,M) 

IMAX  T  =  I DUMR 
IJETT=I JET(M) 

IJETP1=I JET(M)+ 1 
IDUMLL* IOUML 
GO  TO  13 

C*THIS  SECTION  HANDLES  REFRACT/OIFF  IF  THE  T  A0>0 . 0 
12  CONTINUE 

C-FIRST,  REGION  2-  ALL  REFRACTION 
NPTS  =0 

DO  1101 = I DUML . MET (Ml 
J2( I)=J1< I) 

110  J1< I )= 1 

00  111  I  =  I DUML .  I UE  T ( M ) 

DO  111  0  =  0 1 (  I  )  . 02 (  I  )  -  1 

111  NPTS  =  NPTS-M 

IMAXT* I  DUMP 
IDUMLL* IDUML 
IOETT  =  IUET(M1 
I  JE  TP  1  =  I JE  T ( M  )  * 1 

CALL  REFRACI01 . J2 . NPTS ,  IDUMLL . I  JETT . IDUMLL ,M) 

IMAXT  * IDUMR 
IOETT  =  I  JET(ri) 

IOETP 1  *  I  JET ( M 1* 1 
I DUML  L  *  I OUML 

C*NOW  WILL  DO  REGION  3  OF  THE  POS  T HE T AO  CASE 
00  112  I=IJET(M)* 1 . IDUMR 
J2  (  I ) = J  1  (  I  ) 

112  J 1 ( I  )  *  1 

00  113  I=IJET(M1* 1 , IDUMR 
J 1  RE F (  I  )=  1 

C * WI L L  GO  ONE  PT .  BEYOND  J2(I)  TO  MAKE  SURE  OUTOF  DIFF  ZONE 
DO  1 14  J  =  J1( I  )  . J2( I )  +  1 
XCOOR  =  ( 1-10) *DX 
YC00R=0  5*(Y( I , J)+Y( I , J+1 )) 

ANGLE  *  AT  AN ( ( XCOOR -XD I STN )/( SJETTY  YCOOR ) ) 

I F ( YCOOR  GT  SJETTY  )  ANGL E =PI + ANGLE 
C-IF  LEAST  J-VALUE  IS  OUT  OF  SHAD  ZONE , SO  ARE  OTHER  J'S  (FOR  EACH  I) 
I F ( ANGL E . GT  ABS ( THE T AO ( M ) ) )  GO  TO  115 

RAD=SQRT( ( XCOOR -XDI STN)* *2* (SJETTY- YCOOR )**2) 

RHOND  =  RAD*  TWOP I / E  L  T I P 
THE  *  THET  A0( M ) 

CALL  OIFF(RHOND , THE . ANGLE . AMP ) 

ANGRAD=ANGLE 

C*WILL  NOW  REFRACT  DIFFRACTED  WAVES  IN  SHAD  ZONE  USING  SNELL'S 
CTIP=ELTIP/T 

ALPHAS) I ,J)  =  ATAN( (0  5* ( Y ( I  *  1 , J ) ♦ Y ( I  +  1 . J* 1 ) ) -0  5* 
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1 16300 
1 16400 
1 '6500 
1 16600 
1 16700 
1 16800 
1 16900 
1 17000 
1 17100 
1 17200 
1 1 7  300 
1 17400 
1 17500 
1 17600 
1 1 7700 
1 17800 
1 17900 
1 18000 
118100 
1 18200 
1 18300 
1 18400 
1 18500 
1 18600 
1 18700 
t 18800 
1 18900 
1 1 9000 
119100 
1 19200 
1 19300 
1 19400 
1 19500 
1 19600 
1 19700 
1 19800 
1 19900 
120000 
120100 
120200 
120300 
120400 
120500 
120600 
120700 
120800 
1 20900 
121000 
121100 
121200 
121300 
121400 
1 2 1 500 
121600 
1 2 1 700 
1 2 1 800 
121900 
122000 
122 100 
122200 
122300 
122400 
122500 
122600 
122700 
122800 
122900 
123 000 
123 100 
123200 
123300 
123400 
123500 


•  (  Y  I  I  i  ,  .1 1 » i  (  [  i  .  j*  U  1  I  (2  *p»  i  l 

[Fit  EQ  t  JETI  M  )‘  1  IALPHASI  I  .  d  )  ■- A  I  ANI  I  P  5MMI-'.  i  I  •.  M  •  1  1  1 

-  <  v  (  I  . o  )  ♦ yf  I  . i  )  )  )/D*  ) 

DALPHA =ANGRAD  ALPHAS!  I  J) 

THETA!  I  , J ) “ARSINI  ( C 1  I  . J )/CT IP ) ‘SIN! PA!  PH A  1  ! 

THE  I  A!  I  .  0  )  =  THE  1  A!  I  .  J  >  *  ALPHAS!  1  .  .)) 

HI  I , J)=HINC*AMP 

C'MUST  CHECK  TO  SEE  IF  WAVE  WOULD  HAVE  BROKIN 

IF(HBII.J)  LE  HI  I  ,  J  J  AND  HB([,J‘1|  Gf  i«  I  <  I  I  IRRFAK I  I  )  • 

I F I HB (I.JI.LT  HII.J))  HI  I . j I =HR<  I  . HI 
1 14  CONTINUE 
GO  TO  113 
1  15  J 1  REE ( I  )*J 
1 13  CONTINUE 

C • NOW  MUST  DO  REERAC  FOR  RFGION  1 
NPT  5-0 

DO  116  I  =  I JE  T I M ) 4 1  .  I  DUMP 
DO  116  J=  J 1  RE  F I  I  i  . U2 (  I  i  1 
1 16  NPTS=NPTS* 1 
IMAXT= IDUMR 
I DUMLL  =  I DUML 
I  JE  T  T  =  I  JE  T  (  M  1 
I  JE  TP  1  =  I  JE  T  I  M  )  * 1 

CALL  REERACI  J  1RF  F  .  ,)2  .  Ni  I  S  I  ;  MA  .  •  •  m  .<  .  ■  ui 

I  MA  X  T  -  I  DUMR 
I  JE  T  T  =  I  JE  T  I  M  ) 

I JE  T  P  1  -  I  JE  T ( M  1  ‘  1 
IDUMLL= IDUMI 
13  CONTINUE 
200  CONTINUE 
RETURN 
END 

(2*  +  *********»»**»#*««»w»»*»*»»»»». 

SUBROUT  INE  LOCI  I  M.  JJ  .  JO  I  M  .  .JSIM  .  ,  BAR  .  ItlUM  ! 

COMMON/ A/  C( 60. 20) . RKI 60. 20) , Y ( 60. 20 1 . DEEP! 60 . 20  I . ALPHAS!  6.  2  I 

COMMON/ A  A/V7ER01 60) 

COMMON/B/  THETA! 60. 20 ) . Ox TOTI 60) ,  OLDANGI 60 . 20  I .  0< (60.20 
COMMON/C/  H(  60. 20)  .CGI  60. 20  )  .HOLD  I  60.20)  .  HB  I  60. 20  I  .  ■  B  I  FU"  I 
COMMON/N  USED/ JUSE  .  T  .CO  .CC.EN  .  CGGEN  .  ANGC.EN  .  Ox  .  BERM  .  the  t  A0<  <  ■>  I  mma  . 
COMMON/D/ S I GMA . G . ELO . UMAX  .IMAX.PI  TW0PI.PI02. HGEN .  I Jf M  ip  1 . S .. F * ’  - 
C  *  SUBROUT  I NE  I.OC  FINDS  J-VAIUFS  WHICH  ARE  GREATER  AND  IFSS  than  -PAR 
JO  I M  =  2 

2  AA=0  5*(Y( IM. JOIM)+Y( IM, JOIM-  t ) ) 

I F  (  AA  GT  YEAR  |  GO  TO  4 
JO  I M  =  JO  I M+ t 

C  *  THE  FOLLOWING  IS  REDD  SO  THAT  DY/px^O  5 
C*WILL  DTERMINE  K  SIN  THETA  ON  IM-lINE  AT  A  DIST  .BAR 
C'WILL  CALL  THIS  POINT  JUSE‘1 
I F ( JOIM  LE  JUSE I  GO  TO  2 
JOIM  =  JUSE-M 
Y I IM . JOIM ) =YBAR 

C*  DEPTH  AT  THIS  POINT  WILL  BE  COMP  ASSUMING  CONST  BEACH  SLOPE  ON  I  -  I M 
DEL  =  .5* I  V  I  IM . JO IM -  1  ) ‘ Y (  IM . JOIM- ? )  )  5  •  I  V  I  I M . JO  I M -  2 ) *  r I  I M . JO  I M  ill 

BSLOPE  =  (  DE  EP  I  IM .  JO  I M- 2  ) -DEEP!  I M  .  JO  1  M  -  3  )  (/DEL- 

DEEP!  IM.JOIM-  1  )  =DE  EP(  IM  .  JOIM  -  2  )  ‘BSLOPE  •  I  >  I  IM,  JOIM)  v  (  IM  ,  JOIM  .11 
OEPTH=OEEPl IM.JOIM- 1 ) 

C  A  i_  L  WVNUMI  DEPTH,  T  ,  DUMK  ) 

RK ( IM.JOIM- 1 ) =  DUMK 

C( IM, JOIM- 1 )  ~C0*  T  ANH( RK I  IM.JOIM-  1)*DEEP(  IM.JOIM  •)) 

EN  =  0 . 5  *  (  1  0  ‘  (  I  2  0*RK(  IM.JOIM  1  1‘OEErI  IM.JOIM  Ul'MNHI 
*  2  *RK(  IM,  JOIM  -  1  )*E>EEP(  IM.JOIM  -  1  )  )  )  ) 

CGI IM.JOIM- 1 I-CI IM.JOIM  1)*EN 
C  * W I L  L  USE  SNELL'S  LAW  TO  DETERMINE  THE  WAVf  ANGt E  HERE 

C'ANGLE  OF  CONTOUR  WILL  BE  ASSUME  10  BF  IMF  SAMF  AS  IMF  JMAV.i  .  pmi jur 
IFIIOUM  EO  1 (ALPH-ATANI I Y( IM.JOIM  I)  v(IM  1.U0IM  11)0*) 

IFlIOUM  EO.  -  t  )ALPH-ATA‘J((Y(  IM*1.  JOIM  1)  YIIM..IOIM  1)1  n«) 

0 ALPHA  - ANGGEN -ALPH 

THEFAf  IM.JOIM  -  1  )^ARSIN(  (  Cl  IM.JOIM  1  )/rr,FN)  *S  INI  DAI  PHA  1  I 
THETA! IM.JOIM- 1 )*THETA( IM.JOIM  1 )‘Al PH 
4  JS I M  =  JMAX -  1 

6  AA-0.5*(Y<IM,JSIM)‘(Y(IM.JSIM*U)) 

IFIAA.LT. YBAR  )  GO  TO  8 
JSIM= JSIM- 1 


8? 


123600 
123700 
123800 
123900 
1 24000 
124100 
124200 
124300 
124400 
124500 
1 24600 
1 24700 
1 24800 
124900 
125000 
125100 
125200 
125300 
1 25400 
125500 
1 25600 
125700 
125800 
125900 
126000 
126100 
126200 
126300 
126400 
126500 
126600 
126700 
126800 
126900 
127000 
127100 
127200 
127300 
127400 
127500 
127600 
127700 
127800 
127900 
128000 
128100 
128200 
128300 
128400 
128500 
128600 
128700 
128800 
128900 
129000 
129100 
129200 
1 29300 
129400 
129500 
129600 
129700 
129800 
129900 
130000 
130100 
130200 
1 30300 
1 30400 
1 30500 
130600 
130700 
1 30800 


C -  I F  JS 1 M=0 , THERE  IS  NO  ADO  PT .  SUB  REFRAC  CAN  HANDLE  IT 
I F ( JSIM. EO  .  O )  GO  TO  8 

GO  TO  6 
8  RETURN 

END 

C . * . * . 

SUBROUTINE  WVNUMI DEPTH, T . RK ) 

G  =  32  17 
E PS  =  0  OOI 
TW0PI=6  283185307 
S I GMA  =  T  WOP  I /T 

RK  =  T  WOP  I  /  (  T  -  SORT  (  G*  DEPTH  )  I 
DO  100  I T  = 1 , 20 
ARG=RK -DEPTH 

EK  =  ( G-RK*  TANHI ARG )  ) - ( SIGMA  *  *2 ) 
fKPR=G-(ARG-(< SECHI ARG)  I  *  *2 )*TANH( ARG )  ) 

RKNEW'RK-EK/EKPR 

ieiabs(rknew-rk)  le  abs( eps-Rknew  i  >  go  to  120 

RK  =RKNE  W 
*00  CONTINUE 

WRITEI6. 1000)  IT, DEPTH. RK 

1000  FORMAT (///. 10X I TERAT ION  FOR  K  FAILED  TO  CONVERGE  AFTER" 
•  .  3T  .  I  3 .  "  I TfRAT ION"  . / .  "OUTPUT  DEPTH.  RK".3X.2F13  5) 


CALL  EXIT 

1  20 

RK  =RKNEW 

IF ( RK  GT .00) 

GO  TO  140 

WR I TE ( 6 .  1020  ) 

DEPTH . RK 

OUTPUT  DEPTH. RK" . 3X . 2F  1  3  5) 

1020 

FORMAT ( ///  .  10X  , 
CALL  EXIT 

"  RK  IS  NEG" , / 

140 

RETURN 

END 

SUBROUTINE  SMOOTH! THETA . I  MAX , JMAX .  I  JET , S JETTY . MMA  X , Y ) 

C-THIS  WILL  SMOOTH  the  wave  angle  field  to  acct  for  DIFF ( ART  I F ICI ally  ) 

DIMENSION  temp; 60. 20) . Y( 60, 20) . THETA! 60 . 20 ). I  JET (  10 ) 

C  *  f  MMA  X  + 1  )  IS  REQ  D  BECAUSE  M-GROINS  HAVE  M+ 1  REACHES  OF  SHORELINE 
DO  10  M* 1  . MMA X + 1 
IFIM.NE  1 )  GO  TO  3 
I  L  E  F  T  =  2 

IRIGHT=I JET( 1 ) 

GO  TO  5 

3  I F ( M  NE  MMAX+1 )  GO  TO  4 

I  left* i jet ; mmax ) * i 

I R I GHT  =  I MA  X  -  1 
GO  TO  5 

4  ILEFT*IJET(M- 1 )  + 1 
I R I GHT -  I JE  T ( M  ! 

5  CONTINUE 

DO  1  J  = 1  , JMAX -  1 
DO  1  I -ILEFT , IRIGHT 

I F ( I  NE,  ILEFT  AND  I  NE. IRIGHT)  GO  TO  '5 
C*TO  GET  HERE.  MUST  BE  ON  BOUN  OR  ADJ  TO  A  STRUCTURE 
IF ( I . EO . 2 .OR . I . EQ  IMAX- 1 )  GO  TO  15 
C  *  TO  GET  HERE. ADO  TO  A  STRUCT  AND  CAN  BE  ILEFT  OR  IRIGHT 
I F ( Y ( I . J  )  GE  SJETTY)  GO  TO  15 
C  *  I F  HERE.  WITHIN  JETTY  AND  ADJ  TO  EITHER  SIDE. 

IF( I  EO. ILEFT )TEMP( I , J ) =0  5  * ( THE T A ( I . J )  +  THE  T  A ( 1+1 , J) ) 

I F (  I  EQ .  IRIGHT  )  TEMP(  I  , J ) =0  5- < THETA(  I  . J )  +  THET A (  I  -  1  . J )  ) 

GO  TO  1 

15  TEMP ( I . J ) =0  25* THE T  A( I  -  1 , J ) +0  50* THET  A( I . J ) +0  25* THE  TA( I ♦ 1 . J ) 

1  CONTINUE 
10  CONTINUE 

DO  2  J- 1 . JMAX- 1 
DO  2  1=2. IMAX- 1 

2  THETA( I . J)=TEMP( I . J) 

RETURN 

END 

. . . . 

FUNCTION  SECHIA ) 

SECH  = 1  O/COSHIA) 

RETURN 

END 

C****HERE  IS  WHERE  THE  IMSL  ROUTINES  MUST  GO! 
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APPENDIX  C 


CONTOURS  AND  SCHEMATIC  ILLUSTRATIONS 


This  appendix  presents  tables  of  the  original  contours  at  Oregon  Inlet  ai 
the  final  contours  for  the  eight  numerical  simulations  (Tables  C-l  to  C-9) 
Also  included  are  schematic  illustrations  of  sediment  volumes  transported 
from  the  nourished  region  (Figs.  C-l  to  C-8). 


Table  C-l.  Initial  bathymetry  for  all  simulations  (prior  to  any  sediment  addition). 


I m  t*«  CONTOUR  »»iuts.  »,  FUUOa 
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N»Tlrt«  7  7  7  4  a  O  7  0  0  aMM  COO 

9  ;  *  —  3*  r-  *-  7  —  —  «—«**—■—  — 

*»03*r«  —  «r  *o 

CiO»  7  O  —  —  —  O  •«  7  7  -  o  71  7  C  O  •O  i’  "W  O  O  — O  ^  O  O 
•  <V  «  -N  *C  7  m  •»-  3  ^  »NIO|7  »OK  •rf'  O  > 

«|C  OR  -4  O  rf*  »-  K  O  7  «  O  Oi  —  C  -1  O  ^  'M  O  X  *V  —  O  O  *\  7  C 

“ N  •  •  •  «T*  •  «  -  O  «  •  »7  •  •  «C  •••«■»■  •  •  •«>  •  *  «S  •  «  • 

oi-oojo^N.-^-tr-arot  to*>  c  -*  «m- «  mm»»k  ^ 

—  O  O  »  0  9  9  h>  r»  O  ^OO-O/i^-9  o< 

_ ,  9*  9*  n»«T|«V  -1  OO  999  999  999  \\<M  O  C  C 

•i  ^  —  o  t*  —  — - 


U  .*42  «,N.*£  9  9 .5  9  **  *1  ~  - 

♦  *•••?  ••-RIO  499*1  •—•***  . 


M  «C  *  €•  C  -*  O  C  -  « 


«■*  •  •  Rg»  •  •  •  O  •  •  ••  •  •  •«  •  •  »C  •••9  ••  •©  ••• 

J  •  3  «2»o  «»•  *|**5AknO  ^  •  O  —  «  •« 

I  9»  9l9i  AlMN  »999  9  9  9  9*9  9  «l<VM  999 

m  |9  I  9  «•  9  9* - —9  —  —  — 

jl  N  >  •  ,  «4  r  -*  -v 

<  09' 

•  MMM  ••  « 9  •909  9999  •—99'  *99  -  «C  *  (V  *9#*^ 

I  ^  •««*K»N9f  9*9  9(  9*  —  *9  99  0*9140  OR«Of»  9  O  09 
H  •  •  4(0  •  4  *0  4  •  ••!  •  •  #40  •  •  *4D  •  •  •#■►  •  •  •  O  •  •  •» 

^  400*0i*»*  go  ^  9|  0  9*019*  9  9  940  949*0  99  9«|9*»  #•*•  90  «i 

5£jq  jqj  Sj££  |o39  900  Sot  ^iSS^oool 


m 

AI 

« 

9 

— 

0 

0 

0 

p«* 

m 

AI 

— 

0 

O 

nj 

o  >  ^ 

—  n  0 

•  9  **a  • 

—  9 -A- 

mom 

>»  A  O 

0  0  9 

•  4  »a 

•  -•0  . 

•  N  -* 

•  9  9 

•pa  m 

• 

•m  »a 

•  m  ® 

•4  ■«  4) 

n  n  p* 

9  9  9 

v«  n 

O  9  -p» 

0  un  a 

>  *  - 

0mm 

ii  •  • 

AO  •  • 

m  •  • 

o  •  • 

®  ■  ■ 

®i  •  • 

r*  •  m 

9  •  • 

M  -0«« 

AI  O'AI 

**pa<«a 

n9  9- 1 

—  >Am 

-*  *«“  ® 

A»«*.0i 

0  »  n 

0  *a 

9  P* 

P»  9 

m©i 

AI  9 

•♦m  9 

—  m  •• 

—  09 

M  'V 

Ai  AI 

;*A*A 

1  ^  *}• 

*  P* 

ni  — 

A  A 

0  0 

0 

9 

m 

;<m  1 

P» 

•«  —  — 

f—  ••  — 

O  —  — 

Ai 

O 

9. 

1  9  i 

o 

•A 

— 

— 

9  -A  P- 

«  ®  in 

9  9  ni 

m  •  m 

n  mi  o 

(MOtf 

9  *>  m 

m  m  *a 

•1^  *4 

•  —  nj 

•9  « 

•  —  & 

•0  '•« 

■  A® 

•  A  P- 

•  m  — 

•  m  « 

aj  ®  9 

-A  AI© 

le  —  « 

*a  o  pa 

fAj9  0 

m  p*  m 

0  0  9 

9  •  • 

«  •  • 

m|  •  • 

;»  •  • 

r*.  •  • 

mi  •  • 

n<  •  • 

9  •  • 

aj  pa  — 

ns  *>  m 

*\  9  m 

'909. 

n-  91  k» 

—  r-»  ru 

/A  r-  — 

0  om 

•0  9 

9  P- 

If*-  9 

1  PA0 

m  p- 

— m  0 

—  m  ru 

—  09 

nj  -m 

*v  ai 

j  m  *n 

i/i  *' 

®,p^ 

i«j  -« 

fA  pA 

0  0 

m 

<0 

-o 

1  o 

» 

fA  —  — 

•A  —  — 

— *  — 

0 

9 

— 

0 

Ai 

O 

9  i 

— 

0  *a  -. 

o  -a  ®  ao  *a 

<*•4  O  ^ 

'Ao®  —co  in  —  o  ® 

A>  m  o 

0  9  *A  9|  9 

O  'VA 

•  9  in 

o 

•  0  m  m 

•  0  K 

9,  *m  r»  nj 

•  9io  m 

•  PA  ® 

^  »m  ®*  -• 

•  0  9 

-0-0  0 

o 

O  »  7  4  O  —  0 

9i  9  O'  in  o 

0  m  m  r-  9'  O  r\j 

**-  Oi/i  —  9  »*-  om 

<J 

C\J 

O) 
l o 
*0 
U 


to 

=3 

O 


o 

o 


i)  fl  3  •«  o  5.10s  0**9(^9^K4)9 

•  •  •  ■©  •  •  •  •  *N  •  •  *1^1  •  •  tr«l  •  *  *o  •  •  tO  •  •  • 

—  ^  9  OJ  «9>  N  — >  "A  «A  OJ  r\i  n  0  M  p-  *A 

•O'/'-  >  O  O  <\|  —  9  •*OOm^M'0«4|T9^  m  O  to  Q 

'X  ni  m  m  oj  'vj  !  »a  *n  «a!  n  in  9  <a  <a  p*.  I  ru  f\«  o  pa  pa  m  0  0  0 

m  0<  9  Ai  1  r*.|  —  ,M^*cp.««*«w«cr  —  —  — 


,  m  p*  x>  j  pa 

•  m  ®  1  91  pa  ® 

P-^)/(l>il/»iAl«^-'OOi^  co|p«»  Ai  Ni 
I  ••’i  ^  •  —  p-0  «tn  —  f-  •  ®mi  —  •*' 
«  m  0  p*-  —  f*  aj  9  .n  r»  i  9  a  .©i  pa  -•  9(0  n<  ® 


oj  0  0 

9  m  ®|Ai9  -a  9  —  n*  (7-  —  9im  o 
9  «  •  9  9  <7  •  Aim:ni  •«.«  o 
W  9  9\  9  9  0  O  Al  ©Im  ~  O  0 


—m  *a  o  0  Aim  m  —  9  •■ 
r*  ®  0  ««o  —p-  »9ii fs  * 

n#  fM  •  *.  **',*-  m  m  m  M 

I»  •  •  •  •  o  •  »  •■ 

»  «^(V^#OP  A<7  m  • 

r-  m  *a-  •*  »a  mi  -•  •*  0io  •  ■ 

—  Ai  —  »A  »Ai  fA  «,4  4 


88 


Table  C-5.  Final  contours,  case  2.c2. 
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Table  C-6.  Final  contours,  case  2.c3. 
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Table  C-8.  Final  contours,  case  3  (17  weeks  plus  sediment  addition). 
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“as  a  ;<a. 

Period  Considered:  Twelve  months,  Januarv  through  .\cer.ber,  jsir.c  >-;■ 

'.VIS  wave  hindeasts 

Sediment  Budget  Summery: 

.Amount  o:  seoiment  added:  None 

Amount  or  sediment  transported  shoreward  iron  nourished  region:  992  yd' 

Amount  ot  sediment  transported  seaward  iron  nourished  region:  g6  yd- 

Net  amount  or  sediment  transported  alongshore  :ron  nourished  reg :  o vwiq  ( 444  yd' 

local  amount  or  sediment  transported  :  ror.  nourished  region:  -9,356  vd‘ 


Figure  C-l.  Schematic  illustration  of  sediment  volumes  transported 
from  region,  case  2.  a. 
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Figure  C-3.  Schematic  illustration  of  sediment  volumes  transported 
from  nourished  region,  case  2. cl. 
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MICROCOPY  RESOLUTION  TEST  CHART 
NATIONAL  BUREAU  of  STANDARDS- 1965  A 


17-ft  L«v«l 


118,864  yd3 

414.168  yd3 

»  28,920  yd3 

1 

p  -  yV  a,  ^  s  J  J  — ?  1  j 

J 466,260  yd 

-  .  4-ft  l.v.1  - 

336,876  yd 


3 


637,584  yd 


Case  l.cl. 

Period  considered:  Twelve  months,  April  through  March,  using  1975 
WIS  wave  hindcasts. 


S^i****^  Budqat  Suiwnary: 

.Amount  of  sediment  added :  1,452,  000  yd3  (on  7-  and  U-ft  contours) 


Amount  of  sediment  transported  shoreward  from  nourished  region: 

Amount  of  sediment  transported  seaward  from  nourished  region: 

Net  amount  of  sediment  transported  alongshore  from  nourished  region:  421,428  yd3  (29.0pct) 
Total  amount  of  sediment  transported  from  nourished  region:  916,608  yd3  (63.1pct) 


466,260  yd  (32.1pct) 
28,920  yd3  (2.0  pct) 


Figure  C-4.  Schematic  illustration  of  sediment  volumes  transported 
from  nourished  region,  case  2,c2. 


Case  2.c3. 

Period  considered:  Twelve  months,  July  through  June,  using  1975 
WIS  wave  hindcasts. 


Sediment  Budq 


Amount  of  sediment  added:  1,452,000  yd  (on  7-  and  11-ft  contour) 


Amount  of  sediment  transported  shoreward  from  nourished  region:  415,784  yd^  (28. r,  pet) 

Amount  of  sediment  transported  seaward  from  nourished  region:  23,728  yd3  (1.6  pet) 

Net  amount  of  sediment  transported  alonpshore  from  nourished  region330,400  yd^  (22.8pct) 
Total  amount  of  sediment  transported  from  nourished  region:  769,912  yd^  (53.0pct) 


Figure  C-5.  Schematic  illustration  of  sediment  volumes  transported 
from  nourished  region,  case  2.c3. 
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17- ft 

l«v«l  _ 

^  26,452  yd3 

14 -ft 

l«v*l _ - 

_  7  -ft  1«V«1  _  _  , 


_ _ T  395,556  ydj 


Period  considered:  Twelve  months,  October  through  September,  using  1975 
WIS  wave  hindcasts. 

Sfji>jtflt_Bud9tt  Summary: 

Amount  of  sediment  added:  1,452,000  vdJ  (on  7-  and  11-ft  contours). 

Amount  of  sediment  transported  shoreward  form  nourished  region:  39S,SS«  yd3  (27.2  pet) 

Amount  of  sediment  transported  seaward  from  nourished  regron:  26,452  yd3  (2.0  pet) 

Net  amounc  of  sediment  transported  alongshore  from  nourished  region:  346,  144  yd3  (24.0  pcc) 
Total  amount  of  sediment  transported  from  nourished  region:  772,152  yd3  (53.2  pet) 


Figure  C-6.  Schematic  illustration  of  sediment  volumes  transported 
from  nourished  region,  case  2.c4. 


Case  3. 

Period  considered:  Four  months,  January  through  April,  using  1975 
WIS  wave  hindcascs. 

Sediment  ludgt  Summary: 

Amount  of  sediment  added  363,000  yd3  (on  11-  and  14-tt  contours). 
Amount  of  sediment  transported  shoreward  from  nourished  region: 
Amount  of  sediment  transported  seaward  from  nourished  region: 

Net  amount  of  sediment  transported  alongshore  from  nourished  region 
total  amount  of  sediment  transported  from  nourished  region: 


32,164  yd*  (8.9pct) 
4,112  yd3  (l.lpct) 

:  43,708  yd3  (12. 0  pec ) 
79,984  yd3  (22.0  pec) 


Figure  C-7. 


Schematic  illustration  of  sediment  vol 
from  nourished  region,  case  3. 
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392,000  yd3 
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85, 375  yd3 

- 1 - 

_ _ _ _ _ 11-  ft  1ml _ 

119,813  yd3 

96,950  yd3 

135,350  yd3 

115,243  yd3 

139,391  yd3 

134,680  yd3 

- - 

134,373  yd3 
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1 

! 

7-  ft  1ml _ 

_ 1 _ 

•  3 

275,796  yd 


Period  considered:  Twelve  months,  January  through  December,  using  1975 
VfIS  wave  hlndcasts. 


Sediment  Budget  Suwwaryi 


Amount  of  sediment  added:  1,452,000  yd*  (on  7-,  8-,  9-,  and  10-ft  contours). 
Amount  of  sediment  transported  shoreward  from  nourished  region:  275,796  yd* 
Amount  of  sediment  transported  seaward  from  nourished  region:  392,000  yd* 
Net  amount  of  sediment  transported  alongshore  from  nourished  region:  96,679  yd* 
Total  amount  of  sediment  transported  from  nourished  region:  7$4  475  yd* 


(19. Opct) 
(27. C  pet) 
(  6. 7  pet) 
(52.6  pet) 


Figure  C-8. 


Schematic  illustration  of  sediment 
from  nourished  region,  case  4. 
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APPENDIX  D 


METHODOLOGY  ANO  PROGRAM  LISTING  OF  COMPUTER  PROGRAM  WHICH 
CONVERTS  BATHYMETRIC  DATA  INTO  MONOTONICALLY  DECREASING  DEPTH  CONTOURS 


In  order  to  simulate  prototype  shorelines  (and  In  this  case  to  help 
verify  the  numerical  model  via  Channel  Islands  Harbor  data),  the  (x,  y,  z) 
data  points  must  be  transformed  Into  a  form  suitable  for  use  In  the  model 
(l.e.,  bars  can  not  be  present).  First,  the  batfjymetrlc  data  have  to  be  put 
Into  a  form  with  fixed  longshore  and  offshore  spaclngs  (l.e.,  ax  and  ay  equal 
constants).  This  can  be  accomplished  using  one  of  the  many  available  canned 
programs  which  do  the  Interpolation.  The  problem  is  then  one  of  finding  the 
most  suitable  value  of  the  constant.  A,  in  the  equation  h  = 

However,  as  Is  usually  the  case,  the  exact  location  of  the  shoreline  (h  »  0) 
Is  unknown.  In  addition,  one  requires  the  added  constraint  Is  required  that 
the  volumes  of  sediment  (or  conversely,  the  water  above  the  profiles) 
balance.  The  problem  Is  solved  using  LaGrange  Multipliers  and  the  Newton 
Raphson  technique  for  non  linear  equations. 


The  equation  to  be  minimized  Is 
F(A,ydel p  ydel2,  ...  ydelj,^) 


IMAX  IMAX 

1-1  j-l(hmeasi,j  h  Predi,j> 


(D-l ) 


where  A  is  the  scale  parameter  In  the  equilibrium  beach  profile,  ydeU  are  the 
locations  of  the  shoreline  for  the  IMAX  profiles,  h  Is  the  Interpolated 

depth  from  the  survey,  and  h  d  Is  the  depth  predicted  by  the  equation 


(D-2) 


The  constraint  equation  is  as  follows 
g(A,ydel j,  ...  ydel IMAX ^  ^ 


IMAX  (  y  f  «/,  1 

■  &  “U, A,y  * ,deV  d,l 


IMAX 


2  e  ax  A(yf  -  ydel, 
1=1  a  T 


)S/3  a  y 


me  as 


(D-3) 


where  Vpre( j  Is  the  predicted  volume  of  water  above  the  profile  to  the 
reference  datum,  V^as  Is  the  measured  volume  computed  from  the  survey,  ax 
Is  the  longshore  distance  between  onshore-offshore  profiles,  and  yf  is  the 
distance  offshore  to  the  last  point  on  each  of  the  measured  profiles  (It  was 
a  constant  after  the  interpolation  routine  was  used). 
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LaGrange  Multipliers  procedure  says  to  form  the  quantify  F*  as 

*  F*  =  F  -  xg  (D-4) 

take  the  total  differential  of  equation  (0-4) 

dF*  »  dF  -  x  dg  »  dA  ♦  dtydel^  +  ...  afyHiT^  )d(yde1  IMAX}j 

-  x  (s  *  +  3?y3eT7T  d(*deV  +  -  d(*de1  IMAX>) 

(D-5) 

Rearranging 

0  =  dF*  =(k  •  x  "  +(afydei1)  '  x  afydeijij  d(ydeli}  + 

'  (D-6) 

It  is  clear  that  the  terms  in  brackets  in  equation  (D-6)  must  individually 
equal  zero,  however  this  leaves  (IMAX  +  2)  unknown  (udel  i  =  to  IMAX,  A,  and 
x)  and  only  (IMAX  =  1)  Equations.  The  (IMAX  +  2)th  equation  is  taken  as 
equation  (D-3).  The  following  system  of  equation  then  results: 


IMAX  JMAX 
t  Z 
1*1  j  =  l 


[-2(h 


meas 


1 J 


ydel.)2/3)(yi  ydel,)273] 

i  '  i  J  ' 


IMAX  1  C/-7 

-  x  E  i  ax  (y.  -  yde1.)5/  (D-7-1 ) 

i=l  5  '  1 
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0 


dF 

(JTydel 


da  JMAX 

ir  x  ^yde,il=j=i  C2,h"e“u 

*(2/3  A(yi>j-  ydel 1)"1/3 


*(yi  j  -  yde^)273) 

+  x  AX  A  (yf  -  ydel^2^3 
(D-7-2) 


0  ,  dF 

^yde1IMAX} 


.  JMAX  r 

[2(h 


meas 


IMAX.j 


2/3 

‘  A(y IMAX, j"ydel IMAX*  ) 


■1/3 


( 2/3  A(yIMAXj"ydelIMAX*  ^  +  x  ax  A  (yf  -  ydelIMAX) 


2/3 


( D-7- ( I MAX+1 ) ) 


I  MAX 


5/3, 


meas 


=  I  (3/5  ax  A(yf-  yde1Tr/J) 
i  =  l  r  1 


(D-7-(IMAX+2) ) 


Because  Equations  (D-7)  is  a  system  of  nonlinear  equations,  it  can  not 
be  written  in  matrix  form  as  a  CD]  [x]  -  [E]  system  of  equations  (the 
brackets  denote  matrices).  To  solve  the  equations,  a  Newton-Raphson  Iteration 
technique  for  nonlinear  equations  was  used.  This  is  done  by  differentiating 
each  of  the  (IMAX  +  2)  equations  with  respect  to  each  of  the  unknowns,  the 
resulting  equations  are  then  linear  in  terms  of  Aa,  Aydelj,  .  .  . 

AydeliMAX,  ax  .  The  resulting  matrix  is  inverted  to  obtain  the  A(unknown) 
and  the  quantities  are  added  to  the  original  estimates  to  produce  a  better 
estimate.  This  iterative  procedure  is  continued  until  the  changes  become 
acceptably  small.  The  solution  converged  rapidly.  Generally,  the  first  row 
of  the  matrix  to  be  inverted  is  (an  represents  the  kHl  row  and  the  l^h 
column  of  the  matrix). 


all 


IMAX  JMAX  ,,, 

I  Z  2(y.  .  -  ydel.r3 
i=l  j=l  1 


JMAX 

£ 


ydel 


,)  -1/3(h. 


meas  f  J  ‘  “<*1,3  '  *del 


,2/3, 
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JMAX  4 

ll,IMAX+l  *  3  (yIMAX,j 


y^IMAX1  1/3(h 


meas 


IMAXJ 


IMAX  ^  5/3 

al ,IMAX+2  »  l  ii  ax  (y.  -  ydelT)3/J  ] 
1-1  T  1 
The  second  row  of  the  matrix  is  as  follows: 

JMAX 


(0-8) 


a2,l  =  jfj  [3  h  meaSj  j^yl,j  ~  yde11)_x/J  -  ]  A(yx  .-ydel 1)1/J] 


4  -  *  -  -  *-1/3  *  ,  1/3- 


+  xax  (y^  -  ydelj) 


2/3 


JMAX  4 

a«  ij  =  2  [x  A  h 
2,2  j  =  1  meas 


'J'l.j  -^'l>'4/3*  9*2  (*i,J  *  Xdelj)-^3] 

*  sj 


X  (2/3)  AX  A  (yf  -  ydel  ^ 


•1/3 


a2,3  =  0 
a2^ IMAX+1  =  0 

2/3 

a2,IMAX+2  =  Ax  A  (Xf*  Xdeli) 


( D-9 ) 


The  third  row  is  simply  these  elements  repeated  except  that  the  ones  on  the 
right-hand  side  of  the  first  and  last  elements  are  changed  to  twos,  and  the 
83  3  element  Is  similar  to  the  iz,Z  except  the  ones  on  the  right  hand 
side  become  twos.  The  remaining  column  elements  (l.e.,  those  when  the  k  =  1) 
are  zeroes.  This  process  is  continued  to  fill  the  array,  except  for  the  last 
row. 


The  (IMAX+2)t^  row  is  as  follows: 
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IMAX+2, 2  -  -ax  A  lyf  -  yde]^3 
aIMAX+2,  IMAX+1  =  "ax  A  2 /3 

aIMAX+2,  IMAX+2  =  0  (D-iO) 


The  E  matrix  in  the  CD]  Cx]  =  [E]  system  of  equations  is 


- 


IMAX  JMAX 

i*l  j»l  '  2(h  ^Si  J  '  A(y1J  'ydeli )Z/3)(yij 


measj 
IMAX 


a  *x  (yf  -  ydel  .)5/3] 


JMAX 


E?  =  -  I  I  2(h 


j=!  21  *  meas.  .  "  A^i,j  “  ydelj)273) ( (2)  A  (y 


+  x  (ax  A  (yf  -  ydel1)2/3)] 


JMAX 

EIMAX+1  =  "  ^.Z  2(h 


j=l  measIMAX,j‘A(yiMAX^  ‘  yde?IMAX) 


2y 


'((3!  A  <yi>J  .  ydeljl-1/3)  +  X  (ax  A  (y 


IMAX 


■IMAX+2  3  "  E  *  (U)  AX  A(yf  -  ydelr)5/3)  -  V  1 


The  [0]  [x]  ■  [E]  system  of  equations  was  then  solved,  as  explained 
previously,  by  solving  the  x  column  vector  (which  represents  the  changes  In 
the  unknowns,  aA,  aydel]  ...  AydeljMx  ax),  adding  these  changes  to  the 
respective  variables  and  Iterating  unth  a  final  solution  Is  obtained. 

The  computer  program  which  did  these  calculations  for  the  Channel  Island 
Harbor  simulation  follows.  A  user-supplied  matrix  inversion  routine  is 
required  (Line  37,200). 


too  SRESCT  FACE 

200  c*****«********PR0GRAM  CIH/8VALUE 1 

300  FILE  5( KINO "PACK .TITLE* *CIH4  2076 A* ,FILETYPE*7) 

400  FILE  6<KIND-REN0TE) 

SOO  C*THIS  PROGRAM  USES  THE  INTERPOLATED  PROFILES  OF  CIH. 

600  C*IT  FINOS  THE  LOCATION  OF  THE  SHORELINE.  VDEL  AND  THE  BEST 
TOO  C*F I T  LEAST  SOUARES  *8*  VALUE  FOR  H»BY**2/3 

800  C*USES  LAGRANGE  MULTIPLIERS  TO  CONSTRAIN  THE  V0LUMES(S0  THEY  ARE  EQUAL ) 
900  C* THEN  IT  USES  NEWTON-RAPHSON  ITER  FOR  NON-LIN  EOS 
tOOO  DIMENSION  X( 40) 

I  IOO  0 1 ME  NS I ON  WKAR£A(60O).AMATRX(23.23).BMATRX(23, »> 

1200  DIMENSION  Y < 40 . 20 ) . Z( 40 . 20 ) . YOEL ( 40 > , UBEGIN( 40 ) . VDEL1 ( 40) 

1300  DIMENSION  DYTWO( 40. 20) . DYONE < 40. 20) . OYMTWO( 40. 20 > . DYMONE < 40. 20) 

1400  DIMENSION  OYMFOR( 40.  20)  ,OYFOR( 40. 20)  .  YDONE ( 40 . 20)  .  YDMTW0(40.  20) 

1500  DIMENSION  YDMONE ( 40, 20) . YETWO( 40) . YEONE ( 40) . YEMONE ( 40) 

1600  DIMENSION  YEMTWO( 40) . YEMFORl 40 ) . YEF I VE ( 40 ) 

1700  EXP0N*2./3. 

1800  THIRD*0. 3333333333333333 

1900  C*F IRST  READ  IN  THE  PROFILES  FROM  DISKPACX. 

2000  DO  1  I >1,34 

2100  00  1  J»1, IS 

2200  1  REAO( 5 . 100)  X(I),Y(I,U),Z(1.U) 

2300  IOO  FORMAT! 14X , F6 . 0. FS . 0. F5 . 0) 

2400  C'NOW  WE  MUST  GET  A  FIRST  APPROX  FOR  YDEL 

2500  C*WE  WILL  USE  LINEAR  INTERPOLATION  TO  DETERMINE  IT. 

2600  1BEGIN* 1 

2700  IMAX«2 1 

2800  UMAX* IS 

2900  C*CHANGE  PROFILE  TO  SPAN  1  TO  IMAX(IF  ALREADY  DONE, WON'T  HARM  THINGS) 

3000  I  TEMP  1*1 

3100  I TEMP2* IMAX - I BEG IN* 1 

3200  K--1 

3300  DO  777  I*1.ITEMP2 

3400  K*K+ 1 

3500  DO  777  0*1. UMAX 

3600  V(I,U)*Y( IBEGIN+K . U ) 

3700  777  2(1 . U)*Z( IBEGIN+K.U) 

3800  IMAX-ITEMP2 

3900  OX* 100.00 

4000  DO  2  1*1,  IMAX 

4100  DO  3  0*1. UMAX 

4200  IF(Z(I,U). GE .0.0}  GO  TO  3 

4300  C*F IRST  NEG  POINT  ON  THE  PROFILE  IS  SEAWARD  OF  Z-0.0 

4400  C«  WE  MUST  ALSO  REMEMBER  THIS  LOCATION. 

4500  C*IF  Z( I. 1)<0, CHOOSE  ARBITRARY  PT .  ROUTINE  ITERATES  TO  SOLN. 

4600  ZOUM* 1 .0 

4700  IF(U.NE.I)  ZDUM" Z ( I , U- 1 ) 

4800  YOUM*Y( I , U)-50.0 

4900  IF(U.NE.I)  YOUM* Y( I , U“ 1 ) 

5000  OELV*ZOUM/( (ZOUM-Z( I .U) )/( Y( I . J)-YDUM) ) 

5100  VOEL( I )"YOUM*OELY 

5200  UBEGIN( I )*U 

5300  GO  TO  2 

5400  3  CONTINUE 

5500  2  CONTINUE 

5600  C*THC  VALUES  FOR  Z  ARE  NEG  ON  FILE.  MUST  NOW  MAKE  POS . 

5700  C*THE  Z  VALUES  ARE  ALSO  *10. 

5800  00  35  1*1. (MAX 

5900  00  35  U' JBEGIN( I ) . UMAX 

6000  35  Z< I.U)— Z(I,d)/10.0 

6100  C*MUST  INITIALIZE  *B’  SO  WILL  MAKE  A  FIRST  GUESS. 

6200  C*MUST  ALSO  GUESS  LAM80A  ( XL AMB ) 

6300  8*0. 30 

6400  XLAMB* -2 . O 

6500  DO  10  ITER* 1 . IOO 

6600  C*LET'S  CALCULATE  THE  VOL  OF  WATER  ABOVE  THE  PROFILE . VME AS . 

6700  C*tTS  OUR  CONSTRAINT .BUT  SINCE  YDEL  IS  NOT  KNOWN. A  PRIORI. IT  WILL  CHANGE 

6800  VMEAS'O.O 

6900  00  200  1*1. IMAX 

7000  00  200  i)*  JBEG!N(  I ) ,  UMAX 

7100  1F(0.NE . U8EGIN( I ) )  GO  TO  201 


7*00 

7*00 

7400 

7500 
7600 
7700 
7800 
7900 
8000 
•  too 
8200 
8900 
8400 
8500 
8600 
8700 
8600 
8900 
9000 
9100 
9200 
9300 
9400 
9500 
9600 
9700 
9800 
9900 
IOOOO 
10100 
10200 
10300 
10400 
10500 
10600 
10700 
10800 
10900 

I  lOOO 

I I  lOO 
1  1200 
1  1300 
1 1400 
1  1500 
1  1600 
11700 
1  1800 
11900 
12000 
12100 
12200 
12300 
12400 
12500 
12600 
12700 
12800 
12900 
13000 
13100 
13200 
>3300 
13400 
1*900 
1*600 
1*700 
1*600 
1*900 
14000 
14100 
14300 


VMCA$-VNCA$*OX*t( I,d)*(0.5*(Yd.d)*Y( l  .d* 1 )  )-YDEL< I ) ) 

GO  TO  200 

aot  if(j.eq  jmu)  go  to  202 

VMEAS*VMEAS*OX«O.S*<Yd.d*t)-Vd.d-1))*2d.d) 

GO  TO  200 

202  VM£ AS* VME AS*DX  *2(  I  . d) • ( Y(  I  .d)-0.5*(Y(  |  ,  d )♦ Y (  I , d- 1 ) )  ) 

200  CONTINUE 

C* PRIOR  TO  EOS. COMPUTE  AND  STORE  SEVERAL  VALUES  WE  NEED  OVER  AND  OVER 
C'BECAUSE  COMPUTER  CAN'T  RAISE  A  NEG  VALUE  TO  AN  EXPONENT 
C*MUST  PRESERVE  THE  SIGN. 

00  400  I 1-1. IMAX 

DO  401  dd*dBEGlNdI).dMAX 

ARG 1 ■ V ( 1 1 , dd )-YOEL(lI) 

0VS1GN*SIGN( 1 . . ARG 1 ) 

OV*ABS{ V ( Il.dd)-YDEL(II)) 

OYTWO(II ,dd)*DY**EXPON 
OYONE( II ,dd)*OYSIGN*OY«* THIRD 
OVMTWO( I 1 , dd ) *DY •  •  ( -EXPON ) 

OVMONC (II. dd ) *OYS IGN*0 Y • * ( -  THIRD ) 

OYMf 0R( 1 1 , dd)*DY**( -2 . *£XPON) 

DYFOR( 1 1 . dd) *0¥** ( 2 . *EXPON) 

401  CONTINUE 

ARG2- 1400. - VOEL ( II) 

DSIGN*SIGN( 1 . . ARG2 ) 

DYE  *ABS( ARG2 ) 

YETWOIII )-OYE**EXPON 
YEONEdI  )*OSIGN«OY£**THIRO 
VEMONEt II )*OSIGN*OYE**( -THIRD) 

YEMTWOI II )*OYE**( -EXPON ) 

YEMFOR( II) *OYE  * • ( -2 . »EXPON) 

YEFIVE(II) »DSIGN*OYE  *  * ( 5 . ‘THIRD ) 

400  CONTINUE 

C*LET'S  INPUT  THE  FIRST  ROW  OF  THE  MATRIX.  A 
SUM  18*0.0 
00  300  1 1*1. IMAX 
DO  300  dd*0BEGIN( II). dMAX 
300  SUM1B*SUM18*2 . *OYFOR( 1 1 ,  dd) 

AMATRXt 1. 1 )*SUM1B 
SUML AM-00 
00  305  K-1.IMAX 
SUM  IK *0.0 

DO  306  dd-dBEGlN(K) .dMAX 

306  SUM1K«SUM1K»2. *EXPON*OYMONE<K . JJ ) •< 2<K . dd) -2 . *8* 

•  OYTWO(K.dd)) 

SUMLAM-SUMLAM-O  6*0X*YEFIVE(K) 

305  AMATRX( 1 ,K* 1 )*SUM1K 

AMATRX( 1 , 1 MA  X  *  2  )  • SUML  AM 
C*NOW  THE  MIDDLE  ROWS  OF  THE  AMATRX. 

00  410  LR0W*2.IMAX*1 

SUM2B-0.0 

II-LROW-I 

00  4  15  dd*dBEGIN(U  )  .dMAX 

4 15  SUM2B*SUM2B*2 . • EXPON* 2 ( 1 1 ,dd)*DYMONE< 1 1 , dd ) -4 . *EXPON* 

•  8*OVONfdI  .dd) 

AMATRX( LROW .  1 ) *SUM2B*XLAMB*OX*YETWO( I  I ) 

OO  430  11*1, IMAX 
SUM2Y-0.0 

00  429  dd* JBEGINI II), dMAX 

429  $UM2V-SUM2V*2 . *EXPON*THIRO*B*Z< 1 1 , dd ) *OYMFOR( 1 1 , dd )♦ THIRD* EXPON 

•  *2 . •B*8*0YMTW0( 1 1 . dd ) 

IF((II*1).EQ  LROW )  GO  TO  431 
AMATRX(LROW, 1 1 ♦ I ) *0. 0 

GO  TO  430 

431  AMATRX(LR0W. 11*1 )*SUM2Y-XLAMB*EXP0N*DX*B*YEM0NE ( 1 1 ) 

4*0  CONTINUE 

410  AMATRX (LROW, IMAX* 2 ) *DX*B* YE TWO( LROW- 1 ) 

C*NOW  THE  LAST  ROW  OF  THE  MATRIX  A 
SUMF 1*0.0 
DO  490  11-1, IMAX 

490  SUMFB-5UMFS+0.6*DX*VEFIVE( II ) 

AMATRX ( IMAX*2 , 1 ) *SUMFB 


14300 
14400 
14500 
14600 
14700 
14800 
14900 
15000 
15100 
15200 
15300 
15400 
15500 
15600 
15700 
15800 
15900 
16000 
16100 
16200 
16300 
16400 
16500 
16600 
16700 
16800 
16900 
17000 
17  100 
17200 
17300 
17400 
17500 
17600 
17700 
17800 
17900 
18000 
18100 
18200 
18300 


00  493  11*1. IMAX 

453  AMATRX( IMAX*2 . 1 1 ♦ O- -0X*B* VETW0( 1 1 ) 

AMATRX( I MAX *2 , IMAX+2 )*0.0 
C*NOW  MUST  INPUT  THE  BMATRX. 

SUMF 1 A -00 
SUMF 1B*0  0 
00  455  11*1. IMAX 

SUMF 18'SUMF IB+XLAMB  *0 . 6*0X* VEF I VE( 1 1 ) 

00455  dd*dB£GIN( II). UMAX 

455  SUMF lA'SUMF  1A-2 .*(2(11 . dd ) -B*DYTWO( 1 1 . dJ) )»OYTWO( 1 1 , JO) 

BMATRX<  1 .  1 )*-(SUMF 1 A  - SUMF 18) 

00  460  11*1. IMAX 

SUMF 1 1 -0 . 0 

00  462  dd*dB£GlN( I  I ) . JMAX 

462  SUMF 1 1 -SUMF 1 1*2 .*(2(11 , dd ) -B*DYTWO< 1 1 . dd ) ) *EXPON*B*DVMONE( 1 1 . «M) 
SUMF 1 1* SUMF I l*XLAM8*DX*B*Y£TW0( 1 1 ) 

460  BMATRX (11*1. 1)** SUMF 1 1 
SUMV-O.O 

00  465  11*1. IMAX 

465  SUMV-SUMV*0  6*0X*B*YEF IVE ( 1 1 ) 

BMATRX ( I MA  X ♦ 2 , 1 ) ■ - ( SUMV- VMEAS ) 

C*NEXT  LET'S  CALL  THE  MATRIX  INVERSION  ROUTINE  VIA  IMSL 

CALL  LEQT2F ( AMATRX . 1 . I MAX ♦ 2 . 2 3 , BMATRX . 3 . VK ARE A , I ER ) 

C  *  THE  SOLN  IS  RETURNED  IN  THE  VECTOR  BMATRX 
C*F I NALL 7 .  WE  MUST  UPDATE  THE  X  VECTOR  IN  AX-B. 

8*B*BMATRX ( 1.1) 

XLAMB«XLAMB+BMATRX( IMAX  +  2. 1) 

00  470  I  1-1. IMAX 
470  Y0EL(II)*Y0EL(II )+BMATRX ( I !♦ 1 , 1) 

C ‘CHECK  THE  CRITERION  FOR  COMPLETION 
SUMVEC-0.0 
00  475  II-1 . IMAX 
475  SUMVEC*SUMVEC*ABS(BMATRX( 11,1)) 

I F ( SUMVEC . LT . ( 0 . 1  * ( I MAX *2 ) ) )  GOTO  11 

WRITE (6.*/)  B . I  TER . ( I . YOEL ( I ) , I • I , IMAX ) . XLAMB 
lO  CONTINUE 
I  I  CONTINUE 

C*LET ' S  WRITE  IT  ALL  OUT. 

WRITE(6.*/)  ITER.B.(I.YOEUl).I-I.IMAX) 

STOP 

END 
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APPENDIX  E 


USER  DOCUMENTATION  AND  INPUT  AND  OUTPUT  FOR  PROGRAM  VERIFICATION 


The  computer  program  presented  in  Appendix  B  was  run  on  a  Burroughs 
B-7700  computer.  The  B7000/B6000  series  FORTRAN  language  was  designed  so 
several  existing  programs  written  in  FORTRAN  would  be  compatible  with  minimal 
changes.  It  was  designed  to  be  compatible  with  Fortram  IV,  H  level  and  to 
contain  ANSI  X3. 9-1966  Standard  FORTRAN  as  a  subset. 

Line  37,200  of  the  coding  (see  App.  B)  requires  a  subroutine  from  the 
IMSL  subroutine  package,  LEQT1B  and  its  associated  subroutines.  If  the 
user's  computing  center  has  access  to  this  package  of  subroutine  programs 
they  need  only  bind  them  to  the  program  (note:  copyright  laws  prohibited  the 
inclusion  of  the  IMSL  coding).  If  not,  a  substitute  subroutine  must  be  user 
supplied.  It  must  facilitate  the  solution  of  a  banded  storage  mode  matrix. 

The  program  input  will  be  described  here  using  a  card  deck  set-up, 
however,  the  use  of  diskpack  or  magnetic  tape  input  follows  directly.  Lines 
3100,  4100,  5500,  5900,  6800,  7500,  and  12,900  are  read  statements.  The 
cards  used  for  the  simulation  presented  in  this  appendix  are  shown  in  Figure 
E-l.  The  first  card  contains  the  value  of  WDEPTH,  the  depth  of  water  (in 
meters)  to  which  the  input  wave  conditions  are  to  be  transformed  (a  partial 
list  of  variables  used  in  the  program  is  presented  beginning  on  page  A-8  of 
Appendix  A).  The  format  statements  are  obviously  in  the  program  coding. 

The  second  data  input  card  is  read  by  line  4100  where  the  variables 
SJETTY,  BERM,  SFACE ,  and  DIAM  are  required  (length  of  the  structure,  berm 
height,  shore  face  slope,  and  sediment  diameter,  respectively) . 

Lines  5500  reads  MMAX,  the  number  of  structures  to  be  simulated  (as 
set-up  here,  a  maximum  of  10  structures  can  be  modeled,  however,  appropriate 
changes  in  array  dimensions  would  allow  additions  (structures).  Line  5900, 
which  is  in  a  "DO"  loop  reads  the  lesser  I  grid  value  adjacent  to  where  the 
structure  is  desired.  The  number  of  structures,  WAX,  determines  the  number 
of  data  cards  required  here;  3  structures  require  3  cards  with  the  3  I  grid 
locations  (note,  the  present  configuration  of  the  refraction  and  diffraction 
subroutines  requires  evenly  spaced  structures,  however  this  can  be  altered  if 
necessary). 

The  parameter  ADEAN,  which  represents  the  value  of  A  in  the  equilibrium 
profile  used  is  the  next  value  input  (line  6800).  As  mentioned  previously, 
whenever  possible  a  site-specific  value  should  be  used.  The  space-step  and 
time-step  (DX  and  DELT  in  the  coding)  are  input  next  (line  7500). 

The  last  input  values  are  the  wave  data,  HS,  T,  ALPWIS  read  by  line 
12,900.  This  statement  is  in  a  loop  made  by  the  unconditional  GO  TO  statement 
(line  16,400)  and  the  read  statement.  There  is  an  action  specifier  included 
In  the  read  statement  to  transfer  the  program  to  statement  1000,  thereby  stop¬ 
ping  execution  of  the  program  once  all  the  wave  climate  data  have  been  used. 
The  number  of  data  cards  required  for  this  read  statement  is  dictated  by  the 
length  of  the  simulation  and  the  time-step  used. 

The  input  file  and  output  for  program  verification  follow. 
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Figure  E-l.  Card  deck  input  for  program  verification. 
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